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Abstract 


The 15 yr pulsar timing data set collected by the North American Nanohertz Observatory for Gravitational Waves 
(NANOGrav) shows positive evidence for the presence of a low-frequency gravitational-wave (GW) background. 
In this paper, we investigate potential cosmological interpretations of this signal, specifically cosmic inflation, 
scalar-induced GWs, first-order phase transitions, cosmic strings, and domain walls. We find that, with the 
exception of stable cosmic strings of field theory origin, all these models can reproduce the observed signal. When 
compared to the standard interpretation in terms of inspiraling supermassive black hole binaries (SMBHBs), many 
cosmological models seem to provide a better fit resulting in Bayes factors in the range from 10 to 100. However, 
these results strongly depend on modeling assumptions about the cosmic SMBHB population and, at this stage, 
should not be regarded as evidence for new physics. Furthermore, we identify excluded parameter regions where 
the predicted GW signal from cosmological sources significantly exceeds the NANOGrav signal. These parameter 
constraints are independent of the origin of the NANOGrav signal and illustrate how pulsar timing data provide a 
new way to constrain the parameter space of these models. Finally, we search for deterministic signals produced by 
models of ultralight dark matter (ULDM) and dark matter substructures in the Milky Way. We find no evidence for 
either of these signals and thus report updated constraints on these models. In the case of ULDM, these constraints 
outperform torsion balance and atomic clock constraints for ULDM coupled to electrons, muons, or gluons. 


Unified Astronomy Thesaurus concepts: Gravitational waves (678); Cosmology (343); Particle astrophysics (96); 
Gravitational wave sources (677) 
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The standard model (SM) of particle physics currently 
provides our best description of the laws governing the 
universe at subatomic scales. However, it fails to explain 
several observed properties of our universe, such as the origin 
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theories for physics beyond the SM, or BSM theories for short, 
accompanied by a rich experimental program trying to test 
them. The generation of gravitational waves (GWs) is a 
ubiquitous feature of many BSM theories (Maggiore 2000; 
Caprini & Figueroa 2018; Christensen 2019). These GWs form 
a stochastic background and propagate essentially unimpeded 
over cosmic distances to be detected today, whereas electro- 
magnetic radiation does not start free streaming until after 
recombination. Thus, detecting a stochastic GW background 
(GWB) of cosmological origin would offer a unique and direct 
glimpse into the very early universe and herald a new era for 
using GWs to study fundamental physics. 

Cosmological GWBs can be produced by a number of 
particle physics models of the early universe. Notably, cosmic 
inflation generically produces GWs (Guzzetti et al. 2016), 
which may be observable at nanohertz frequencies if their 
energy density spectrum is sufficiently blue-tilted. Similarly, an 
enhanced spectrum of short-wavelength scalar perturbations 
produced during inflation can source so-called scalar-induced 
GWSs(SIGWs; Doménech 2021; Yuan & Huang 2021а). 
Another potential source of GWs are cosmological first-order 
phase transitions (Caprini et al. 2016, 2020; Hindmarsh et al. 
2021), which proceed through bubble nucleation; bubble 
collisions and bubble interactions with the primordial plasma 
giving rise to sound waves contribute to GW production. 
Finally, topological defects left behind by cosmological phase 
transitions, such as cosmic strings and domain 
walls (Vilenkin 1985; Hindmarsh & Kibble 1995; Sai- 
kawa 2017), can radiate GWs and hence contribute to 
the GWB. 

The North American Nanohertz Observatory for Gravita- 
tional Waves (NANOGrav; McLaughlin 2013) has recently 
found the first convincing evidence for a stochastic GWB 
signal, as detailed in Agazie et al. (20236, hereafter NG15gwb). 
Analyzing 15 yr of pulsar timing observations, NANOGrav has 
detected a red-noise process whose spectral properties are 
common among all pulsars and that is spatially correlated 
among pulsar pairs in a manner consistent with an isotropic 
GWB. In the following, we will refer to this observation as “the 
NANOGrav signal," “the GWB signal,” or simply “the signal,” 
keeping in mind the level of statistical significance at which the 
GW nature of the signal has been demonstrated in NGI5gwb. 
While the GWB is primarily expected to arise from a 
population of inspiraling supermassive black hole binaries 
(SMBHBs; Rajagopal & Romani 1995; Jaffe & Backer 2003; 
Wyithe & Loeb 2003; Sesana et al. 2004; Burke-Spolaor et al. 
2019), cosmological sources may also contribute to it. 

The SMBHB interpretation of the signal is considered in 
Agazie et al. (2023c, hereafter NG15smbh). In this paper, we 
analyze the NANOGrav 15 yr data set (Agazie et al. 2023a, 
hereafter NG15) to investigate the possibility that the observed 
signal is cosmological in nature or that it arises from a 
combination of SMBHBs and a cosmological source. In 
particular, we consider phenomenological models of cosmic 
inflation, SIGWs, first-order phase transitions, cosmic strings 
(stable, metastable, and superstrings), and domain walls. We 
find that all of these models, except for stable cosmic strings of 
field theory origin, are consistent with the observed GWB 
signal. Many models provide in fact a better fit of the 
NANOGrav data than the baseline SMBHB model, which is 
reflected in the outcome of a comprehensive Bayesian model 
comparison analysis that we perform: several new-physics 
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models result in Bayes factors between 10 and 100. We also 
consider composite models where the GWB spectrum receives 
contributions from new physics and SMBHBs. Comparing 
these composite models to the SMBHB reference model leads 
to comparable results, again with many Bayes factors falling 
into the range from 10 to 100. Cosmic superstrings, as 
predicted by string theory, are among the models that provide a 
good fit of the data, while stable cosmic strings of field theory 
origin only result in Bayes factors in the range from 0.1 to 1. 

The reason that some of the Bayes factors reach large values 
is that the SMBHB signal expected from the theoretical model 
used in this analysis agrees somewhat poorly (only at the level 
of 95% regions) with the observed data, leaving room for 
improvement by adding additional sources or better noise 
modeling. It is perhaps an intriguing idea that this disagreement 
may point to the presence of a cosmological source, but the 
present evidence is quite weak. We stress that Bayes factors for 
additional models beyond the SMBHB interpretation are highly 
dependent on the range of priors with which these models are 
introduced. Thus, one should not assign too much meaning to 
the exact numerical values of the Bayes factors reported in 
this work. 

In many models, there are ranges of parameter values that 
would produce signals in conflict with the NG15 data. In those 
cases, we show the excluded regions and give numerical upper 
limits for individual parameters. We do so in terms of a new 
statistical test, introducing what we call the K ratio. These 
parameter constraints are independent of the origin of the signal 
in the №015 data and a testament to the constraining power of 
PTA data in the search for new physics. In our parameter plots, 
we label the K-ratio constraints by NGI5, and where 
applicable, we compare them to other existing bounds. In 
many cases, the №015 bounds are complementary to existing 
bounds, highlighting the fact that new-physics searches at the 
PTA frontier venture into previously unexplored regions of 
parameter space. 

Aside from cosmological GWBs, signals of new physics can 
appear in GW detectors in a deterministic manner. Although 
pulsar timing arrays (PTAs) are primarily used to search for a 
GWB, we can also leverage their remarkable sensitivity to 
search for these deterministic signals. Specifically, DM 
substructures within the Milky Way can produce a Doppler 
effect by accelerating Earth or a pulsar (Seto & Cooray 2007), 
ога Shapiro delay of the photons' arrival times by perturbing 
the metric along the photon geodesic (Siegel et al. 2007). PTAs 
can also probe models of ultralight DM (ULDM), which can 
cause shifts in the observed pulse timing via metric fluctuations 
(Khmelnitsky & Rubakov 2014; Porayko & Postnov 2014) or 
via couplings between ULDM and SM particles (Graham et al. 
2016; Kaplan et al. 2022). We search for both of these 
deterministic signals, and after finding no evidence for either of 
them, we derive new bounds on both these models. 

This paper is organized as follows. We describe the NG15 
data set in Section 2 and our general analysis methods in 
Section 3. In Section 4, we discuss the GWB expected from 
SMBHBs. We present the analysis and results for new-physics 
models that generate a cosmological GWB in Section 5 and for 
models that produce deterministic signals in Section 6. We 
conclude in Section 7. Additionally, we include a list of 
parameters for each model, the prior ranges we use in our 
analysis, and the corresponding recovered posterior ranges in 
Appendix A. We present median GW spectra for all 
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cosmological models based on our recovered posterior 
distributions in Appendix B, and we provide supplementary 
material for specific models in Appendix C. 


2. PTA Data 


The NANOGrav 15 yr (NGI5) data set consists of observa- 
tions of 68 millisecond pulsars made between 2004 July and 
2020 August. This updated data set adds 21 pulsars and 3 yr of 
observations to the previous 12.5 yr data set (Alam et al. 2021). 
One pulsar, J0614—3329, was observed for less than 3 yr, 
which is why it is not included in our analysis. The remaining 
67 pulsars were all observed for more than 3 yr with an 
approximate cadence of 1 month (with the exception of six 
pulsars that were observed weekly as part of a high-cadence 
campaign, which started in 2013 at the Green Bank Telescope 
and in 2015 at the Arecibo Observatory). 

The pulse times of arrival (TOAs) were generated from the 
raw data following the procedure discussed in Arzoumanian 
et al. (2015, 2018а) and Alam et al. (2021). The resulting 
cleaned TOAs were fit to a timing model that accounts for the 
pulsar's period and spin period derivative, sky location, proper 
motion, and parallax. For pulsars in a binary system, we 
included in the timing model five Keplerian binary parameters 
and an additional non-Keplerian parameter if they improved the 
fit as determined by an F-test. Pulse dispersion was modeled as 
a piecewise constant with the inclusion of DMX 
parameters (Arzoumanian et al. 2015; Jones et al. 2017). The 
timing model fits were performed using the TT(BIPM2019) 
timescale and the JPL Solar System Ephemeris model 
DE440 (Park et al. 2021). Additional detail about the data set 
and the processing of the TOAs can be found in NGI5 and 
Agazie et al. (2023d, hereafter NG15detchar). 


3. Data Analysis Methods 


The statistical tools needed to describe noise sources, GWBs, 
and deterministic signals in pulsar timing data have already 
been extensively discussed in the literature (see, e.g. 
Arzoumanian et al. 2016, 2018b). In the following brief 
overview, we focus on the implementation of new-physics 
signals within this framework. 


3.1. Likelihood 


Our search for a new-physics signal utilizes the pulsars' 
timing residuals, ôt. These timing residuals measure the 
discrepancy between the observed TOAs and the ones predicted 
by the pulsar timing model described in NGI5 and briefly 
summarized in Section 2. There are three main contributions to 
these timing residuals: white noise, time-correlated stochastic 
processes (also known as red noise), and small errors in the fit to 
the timing-ephemeris parameters (Vallisneri et al. 2020). Speci- 
fically, we can model the timing residuals as 


ót —n + Fa + M є. (1) 


In the remainder of this section, we will define and discuss each 
of these three terms and define the РТА likelihood. 

The first term on the right-hand side of Equation (1), m, 
describes the white noise that is assumed to be left in each of 
Фе № од timing residuals after subtracting all known 
systematics. White noise is assumed to be a zero mean normal 
random variable, fully characterized by its covariance. For the 
receiver/back-end combination J, the white-noise covariance 
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(nin) = Filo? зум + Qi бу + 7170, (2) 


where i and j index the TOAs, c; s/w is the TOA uncertainty for 
the ith observation, 7; is the Extra FACtor (EFAC) parameter, 
О; is the Extra QUADrature (EQUAD) parameter, and Jr is 
the ECORR parameter. ECORR is modeled using a block 
diagonal matrix, М, with values of 1 for TOAs from the same 
observing epoch and zeros for all other entries. Following the 
approach of previous works (Arzoumanian et al. 2016, 2018b), 
we fix all white-noise parameters to their values at the maxima 
in the posterior probability distributions recovered from single 
pulsar noise studies in order to increase computational 
efficiency (NG15detchar). 

Time-correlated stochastic processes, like pulsar-intrinsic red 
noise and GWB signals, are modeled using a Fourier basis of 
frequencies i/T ps, where i indexes the harmonics of the basis 
and To»s is the timing baseline, extending from the first to the 
last recorded TOA in the full PTA data set. Since we are 
generally interested in processes that exhibit long-timescale 
correlations, the expansion is truncated after М, frequency bins. 
In this paper, we use Ny = 30 for pulsar-intrinsic red noise and 
Мр = 14 for GWBs. The latter choice stems from the 
observation that most of the evidence for a GWB comes from 
the first 14 frequency bins. More specifically, fitting a 
common-spectrum uncorrelated red-noise process with a 
broken power-law spectral shape to the NGIS data, the 
posterior distribution for the break frequency reaches it 
maximum around the 14th frequency bin (NGI5gwb). This 
set of 2N; sine-cosine pairs evaluated at the different 
observation times is contained in the Fourier design matrix, 
F. The Fourier coefficients of this expansion, a, are assumed to 
be normally distributed random variables with zero mean and 
covariance matrix, (aa!) = ф, given by 


llano = bj Tard; + bab Pai)» (3) 


where a and b index the pulsars, i and j index the frequency 
harmonics, and Г,ь is the GWB overlap reduction function, 
which describes average correlations between pulsars a and b 
as a function of their angular separation in the sky. For an 
isotropic and unpolarized GWB, Гь is given by the Hellings & 
Downs correlation (Hellings & Downs 1983), also known as 
“quadrupolar” or “HD” correlation. 

The first term on the right-hand side of Equation (3) 
parameterizes the contribution to the timing residuals induced 
by a GWB in terms of the model-dependent coefficients @;. In 
this work, we consider two kinds of GWB sources: one of 
astrophysical origin, namely a population of inspiraling 
SMBHBs (discussed in Section 4), and one of cosmological 
origin, induced by one of the exotic new-physics models under 
consideration (discussed in Section 5). The last term in 
Equation (3) models pulsar-intrinsic red noise in terms of the 
coefficients р, і, where 


д2 ү г ү, 
EL p 4 
#0? 127? zh 4 "n 


and pai = Yali/Tovs) for all М, frequencies. The priors for the 
red-noise parameters are reported in Table 2. 

Finally, deviations from the initial best-fit values of the m 
timing-ephemeris parameters are accounted for by the term M e. 
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The design matrix, M, is ап Nyo4 x т matrix containing the 
partial derivatives of the TOAs with respect to each timing- 
ephemeris parameter (evaluated at the initial best-fit value), and 
€ is a vector containing the linear offset from these best-fit 
parameters. 

Since in this analysis we are not interested in the specific 
realization of the noise but only in its statistical properties, we 
can analytically marginalize over all the possible noise 
realizations (i.e., integrate over all the possible values of a 
and e). This leaves us with a marginalized likelihood that 
depends only on the (unknown) parameters describing the red- 
noise covariance matrix (1.е., Aa, Ya, plus any other parameters 
describing ; уап Haasteren & Levin 2012; Lentati et al. 
2013): 


1 =T => 
exp (- 5 ot ccs) 
JdetQ2zC) 


where С = № + ТВТ". Here М is the covariance matrix of white 
noise, Т = [М, Е], and B=diag(co, Фф), where оо is a 
diagonal matrix of infinities, which effectively means that we 
assume flat priors for the parameters in e. Since in our 
calculations we always deal with the inverse of B, all these 
infinities reduce to zeros. 

Equation (5) can be easily generalized to take into account 
deterministic signals (like the ones that will be discussed in 
Sections 6.1 and 6.2). In the presence of a deterministic signal, 
h(0), which depends on a set of parameters 6, we just need to 
shift the residuals, ót — ót — h(0). 

Finally, we relate our characterization of the GWB given in 
Equation (3) in terms of Ф; to the commonly adopted spectral 
representation in terms of the GWB energy density per 
logarithmic frequency interval, dpgy/d Inf, as a fraction of 
the closure density, 1.е., the total energy density of our 
universe, pc (Allen & Romano 1999), 


| dogw(f) _ 81? ФСР) 
аў H Af 


p(ot\d) = (5) 


Qew(f) = (6) 


Here Но is the present-day value of the Hubble rate, 
Af=1/Tops is the separation between the М; frequency bins, 
and ®(f) determines the coefficients Ф, in Equation (3), i.e., 
Ф; = Ф(1/Тьѕ). Note that ®(f) is identical to the timing 
residual power spectral density (PSD), S(f) = ®(f)/Af, up to 
the constant factor of 1/Af. In the remainder of this paper, we 
will often work with /^ Ow instead of Осуу, where h is the 
dimensionless Hubble constant, Hp = п x 100 km . Mpc |," 
such that the explicit value of Но cancels in the product ie ie: 


3.2. Bayesian Analysis 


The goal of this work is to investigate a series of 
cosmological interpretations of the GWB signal in our data. 
Specifically, we would like to answer two questions. First, what 
is the region in the parameter space of the new-physics models 
that could produce the observed GWB? And second, is there 
any preference between the astrophysical and cosmological 
interpretations of the signal? 

To answer these questions, we make use of Bayesian 
inference. Bayesian inference is a statistical method in which 
Bayes’s rule of conditional probabilities is used to update one’s 
knowledge as observations are acquired. Given a model H, a 
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set of parameters O, and data D, we can use Bayes's rule to 
write 
P(DJO, H)P(OlH) 


P(O[D, H) = РФ 


(7) 


where P(O[D, H) is the posterior probability distribution for 
the model parameters, P(D|O, H) is the likelihood, P(O[71) is 
the prior probability distribution, and 


ZzP(D[A)- fae P(D|O, H)P(O|H) (8) 


is the marginalized likelihood, or evidence. In the context of 
this work, H is the timing residual model given in Equation (1), 
© contains the parameter describing the covariance matrix @, 
and the data are the timing residuals ôt. The likelihood function 
for our analysis is given by Equation (5) and implemented 
using the ENTEPRISE EXTENSIONS (Hllis et al. 2019) and 
ENTEPRISE EXTENSIONS (Taylor et al. 2021) packages. 
Our prior choices are summarized in Tables 2 and 3. 

The posterior distribution on the left-hand side of 
Equation (7) is the central result of the Bayesian analysis and 
contains all the information needed to answer our two original 
questions. Indeed, integrating over all the model parameters 
except one (two) allows us to derive marginalized distributions 
that can be used to obtain 1D (2D) credible intervals. At the 
same time, given two models Họ and Ні, we can perform 
model selection by calculating the Bayes factor defined as 

21 РИФІНО) 


D) = — = Я 
BoD = = PO) 


(9) 


The numerical value of the Bayes factor for a given model 
comparison can then be interpreted as evidence against or in 
favor of model hypothesis Hı according to the Jeffreys 
scale (Jeffreys 1961): Від < 1 means that Hj is disfavored, 
while By) values in the ranges [10?9, 10991, [10°°, 1019), 
[10 9, 10121, [10'°, 107°], [10?9, оо) are interpreted as 
negligibly small, substantial, strong, very strong, and decisive 
evidence in favor of 7, respectively. 

Given the large number of parameters, the integration 
required to derive marginalized distributions and Bayes factors 
needs to be performed through Monte Carlo sampling. 
Specifically, we use the Markov Chain Monte Carlo (MCMC) 
tools implemented in the PTMCMCSampler package (Ellis & 
van Haasteren 2017) to sample from the posterior distributions. 
The marginalized posterior densities shown in our plots are 
then derived by applying kernel density estimates to the 
MCMC samples via the methods implemented in the 
GetDist package (Lewis 2019). 

In order to compute the Bayes factor between two models, 
we use product space methods (Carlin & Chib 1995; God- 
sill 2001; Hee et al. 2015), instead of calculating the evidence 
Z for each model separately. This procedure recasts model 
selection as a parameter estimation problem, introducing a 
model indexing variable that is sampled along with the 
parameters of the competing models and controls which model 
likelihood is active at each MCMC iteration. The ratio of 
samples spent in each bin of the model indexing variable 
returns the posterior odds ratio between models, which 
coincides with the Bayes factor for equal model priors, 
P(H)- Р(Но). The Monte Carlo sampling uncertainties 
associated with this derivation of the Bayes factors can be 
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estimated through statistical bootstrapping (Efron & Tibshir- 
ani 1986). Bootstrapping creates new sets of Monte Carlo 
draws by resampling (with replacement) the original set of 
draws. These sets of draws act as independent realizations of 
the sampling procedure and allow us to obtain a distribution for 
the Bayes factors from which we derive point values and 
uncertainties on our Bayes factors corresponding to mean and 
standard deviation. Specifically, the central values and 
corresponding errors quoted in the following for the Bayes 
factors were derived by creating 5 x 10^ realizations of our 
Monte Carlo draws. 

From Equation (8), it is evident that models’ evidence and, 
therefore, Bayes factors depend on the prior choice. In our 
analysis, we will often restrict priors to the region of parameter 
space for which cosmological models produce an observable 
signal in the PTA frequency band. However, a more 
appropriate prior choice would cover the entire allowed region 
of parameter space. Nonetheless, when working with flat 
priors, it is easy to rescale the Bayes factors to account for 
wider prior ranges. Specifically, if the priors are extended to a 
region of parameter space for which the likelihood P(D|O, H) 
is approximately zero, the Bayes factors decrease by a factor 
proportional to the increase in prior volume. 

For each model H considered in our analysis, we use the 
reconstructed posterior distribution, P(O[D, H), to identify 
relevant parameter ranges and set upper limits. Specifically, we 
identify 68% (95%) Bayesian credible intervals (Bernardo & 
Smith 2000) by integrating the posterior over the regions of 
highest density until the integral covers 6896 (9596) of the 
posterior probability. Moreover, we give upper limits above 
which the additional model is “strongly disfavored” according 
to the Jeffreys scale (Jeffreys 1961). For instance, to place a 
bound on a single parameter 0, we first marginalize over all 
other model parameters and then determine the parameter value 
at which the likelihood ratio 


gig = 70. (10) 
P(Dl0o, H) 


1 
has dropped to К = 30 Here бо refers to the parameter limit in 


which the new-physics contribution to the total signal becomes 
negligible and P (D|0, H) no longer depends on the exact value 
of 0. Graphing P(D|0, H) as a function of 0, this parameter 
region appears as a plateau, with P(D|@, H) denoting the 
height of this plateau. Assuming a flat prior оп 0, the ratio in 
Equation (10) is identical to the corresponding ratio of 
marginalized posteriors. Furthermore, multiplying and dividing 
by the prior on 6, 


P(0[H) POD, H) 


K(8) = 
P(@I|D, H) PUOH) 


(11) 


The first factor is the Savage-Dickey density ratio and can 
hence be identified as the Bayes factor B = P(D|H)/P(D|Ho), 
where Họ is the model that results from model H when 
omitting the signal contribution controlled by the parameter 0. 
The K ratio can thus be written as the product of the global 
Bayes factor and the local posterior-to-prior ratio for the 
parameter 6, 


P(0[D, H) 


К(0) = В 
m P(0[H)) 


(12) 
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Once В is known, it is straightforward to evaluate 
Equation (12) and determine the K-ratio bound on б. 
Equation (12) is useful for numerically evaluating K, as it 
automatically encodes the height of the plateau in the 
marginalized posterior, P(@|D, H) = P(0|4)/B, which we 
would otherwise have to obtain from a fit to our MCMC data. 
However, we stress that K is defined as a likelihood ratio, 
which renders it immune to prior effects (prior choice, range, 
etc.; Azzalini 1996). For more than one parameter dimension, 
we proceed analogously and derive bounds based on the 


criterion К (©) > us 


All Bayesian inference analyses discussed in this work were 
implemented into ENTEPRISE EXTENSIONS via a newly 
developed wrapper that we call PTArcade (Mitridate 2023; 
Mitridate et al. 2023, in preparation). This wrapper is intended 
to allow easy implementation of new-physics searches in PTA 
data. We make this wrapper publicly available at doi:10.5281 / 
zenodo.7876429. Similarly, all MCMC chains analyzed in this 
work can be downloaded at doi:10.5281 /zenodo.8010909. 


4. GWB Signal from SMBHBs 


Most galaxies are expected to host a supermassive black hole 
(SMBH) at their center (Kormendy & Ho 2013; Akiyama et al. 
2019). During the hierarchical merging of galaxies taking place 
in the course of structure formation (White & Rees 1978), these 
black holes are expected to sink to the center of the merger 
remnants, eventually forming binary systems (Begelman et al. 
1980). The gravitational radiation emitted by this population of 
inspiraling SMBHBs forms a GWB in the PTA band (Rajagopal 
& Romani 1995; Jaffe & Backer 2003; Wyithe & Loeb 2003) 
and is a natural candidate for the source of the signal observed in 
our data. 

The shape and normalization of this GWB depend on the 
properties of the SMBHB population and on its dynamical 
evolution (Enoki & Nagashima 2007; Sesana et al. 2008; 
Kocsis & Sesana 2011; Kelley et al. 2017). As discussed 
in NG15smbh, the normalization is primarily controlled by the 
typical masses and abundance of SMBHBs, while the shape of 
the spectrum is determined by subparsec-scale binary evol- 
ution, which is currently unconstrained by observations. For a 
population of binaries whose orbital evolution is driven purely 
by GW emission, the resulting timing residual PSD is a power 
law with a spectral index (defined below in Equation (13)) of 
-"ївнв = —13/3 (Phinney 2001), produced by the increasing 
rate of inspiral and decreasing number of binaries emitting over 
each frequency interval. However, as GW emission alone is 
typically insufficient to merge SMBHBs within a Hubble time, 
the number of binaries emitting in the PTA band depends on 
interactions between binaries and their local galactic environ- 
ment to extract orbital energy and drive systems toward 
merger (Begelman et al. 1980). If these environmental effects 
extend into the РТА band, or if binary orbits are substantially 
eccentric, Шеп the GWB spectrum can flatten at low 
frequencies (typically expected at f<lyr';Kocsis & 
Sesana 2011). At high frequencies, once the expected number 
of binaries dominating the GWB approaches unity, the 
spectrum  steepens below 13/3 (typically expected at 
f>1 yr Sesana et al. 2008). 

Unfortunately, current observations and numerical simula- 
tions provide only weak constraints on the spectral amplitude 
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or the specific locations and strengths of power-law deviations. 
Despite these uncertainties, the sensitivity range of PTAs is 
sufficiently narrowband that it is reasonable, to first approx- 
imation, to model the signal by a power law in this frequency 
range: 


A&m 1 (_f | "s 
Ф = г, 13 
BHB(S) 127? Ta. | yr! y (13) 


where Фвнв/ Af is the timing residual PSD (see Equation (6)). 

Following Middleton et al. (2021), we can gain some insight 
into the allowed range of values for the amplitude, Авнв, and 
slope, увнв, of this power law by simulating a large number of 
SMBHB populations covering the entire range of allowed 
astrophysical parameters. Specifically, we consider the 
SMBHB populations contained in the GWOnly-Ext library 
generated as part of ће NG15smbh analysis (and discussed in 
additional detail there). This library was constructed with the 
holodeck package (L. Z. Kelley et al. 2023, in preparation) 
using semianalytic models of SMBHB mergers. These models 
use simple, parameterized forms of galaxy stellar mass 
functions, pair fractions, merger rates and SMBH- 
mass versus galaxy-mass relations to produce binary popula- 
tions and derived GWB spectra. While some parameters in 
these models are fairly well known (e.g., concerning the galaxy 
stellar mass function), others are almost entirely unconstrained 
—particularly those governing the dynamical evolution of 
SMBHBs on subparsec scales (Begelman et al. 1980). The 
GWOnly-Ext library assumes purely GW-driven binary 
evolution and uses relatively narrow distributions of model 
parameters based on literature constraints from galaxy-merger 
observations (e.g., Tomczak et al. 2014) in addition to more 
detailed numerical studies of SMBHB evolution (e.g., 
Sesana 2013). 

For each population contained in the GWOn1y-Ext library, 
we perform a power-law fit of the corresponding GWB 
spectrum across the first 14 frequency bins that we use in our 
analysis. The distribution for Авнв and увнв obtained in this 
way is reported in Figure 1 (blue contours) and compared to the 
results of a simple power-law fit to the GWB signal in the 
NGI5 data set (green contours). The 95% regions of the two 
distributions barely overlap, signaling a mild tension between 
the astrophysical prediction and the reconstructed spectral 
shape of the GWB. In view of this observation, we stress again 
that while these simulated populations are consistent with 
systematic investigations of the GWB spectrum (e.g., 
Sesana 2013), they assume circular orbits and GW-only driven 
evolution. Adopting models that include either significant 
coupling between binaries and their local environments or very 
high eccentricities could serve to flatten the spectral shape and 
lead to SMBHB signals that better align with the observed data 
(see NG15smbh for an extended discussion). Neither of these 
effects, however, is expected to significantly impact the 
amplitudes of the predicted spectra that, for expected values 
of astrophysical parameters, remain in mild tension with 
observed data. As discussed in NGl5smbh, in order to 
reproduce the observed amplitude, SMBHB models require 
one or more of the astrophysical parameters describing the 
binaries’ population to differ from expected values. For the 
present analysis, the spectra derived from the GWOnly-Ext 
library thus represent a convenient benchmark that is simple, 
well defined, and easy to use. By using theory-motivated 
priors, our reference model constitutes an important step 
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Figure 1. Comparison of the 68% and 95% probability regions for the 
amplitude and slope of a power-law fit to the observed GWB signal (green 
contours) and predicted for purely GW-driven SMBHB populations with 
circular orbits (blue contours; МС 1 55тЬћ). The black dashed lines represent a 
2D Gaussian fit of the blue contours. The vertical red line indicates y = 13/3, 
the naive expectation for a GWB produced by a GW-driven SMBHB 
population (Phinney 2001). 


toward a more realistic modeling of the GWB spectrum from 
inspiraling SMBHBs that goes beyond a power-law parameter- 
ization with spectral index "внв = 13/3, which has been the 
standard reference model in much of the PTA literature over 
the past decades. 

The black dashed contours in Figure 1 show the results of a 
2D Gaussian fit to the distribution of Авнв and ygyp values 
derived from the simulated SMBHB populations (see 
Equation (A1) in Appendix A for the parameters of this 
Gaussian distribution). This fitted distribution is what we adopt 
as a prior distribution for Авнв and увнв in all parts of the 
analysis described in this paper. 


5. GWB Signals from New Physics 


In this section, we discuss the GWB produced by various 
new-physics models and investigate each model alone and in 
combination with the SMBHB signal as a possible explanation 
of the observed GWB signal. For each model, we give a brief 
review of the mechanism behind the GWB production and 
discuss the parameterization of its signal prediction. We report 
the reconstructed posterior distributions of the model para- 
meters and compute the Bayes factors against the baseline 
SMBHB interpretation. In Figure 2, we show a summary of 
these Bayes factors; in Figure 3, we present median 
reconstructed GWB spectra in the PTA band for a number of 
select new-physics models; and in Figure 4, we show similar 
median reconstructed GWB spectra in the broader landscape of 
present and future GW experiments. 

As discussed in Section 4 and in more detail in NG15smbh, 
there is a mild tension between the NGI5 data and the 
predictions of SMBHB models. The models generally prefer a 
weaker and less blue-tilted hi ew spectrum than the data. This 
discrepancy presents an opportunity for new-physics models to 
fit the data better than the conventional SMBHB signal. 
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Eventually, this tension may grow to the point of giving strong 
evidence for new physics, or it may be resolved with better 
modeling and more data. Specifically, models of SMBHB 
evolution with a significant coupling between binaries and their 
local environment could lead to a signal that better aligns with 
the data and reduce the evidence for new physics. For all these 
reasons, we caution against overinterpreting the observed 
evidence in favor of some of the new-physics models discussed 
in the following sections. 


5.1. Cosmic Inflation 
5.1.1. Model Description 


Cosmic inflation denotes a stage of exponential expansion in 
the early universe that provides an explanation for the initial 
conditions of big bang cosmology (Liddle & Lyth 2009). At the 
level of the background expansion, inflation accounts for the 
size, homogeneity, isotropy, and flatness of the observable 
universe on cosmological scales; at the level of perturbations, it 
provides the seeds for structure formation in the form of 
primordial density fluctuations. In the standard scenario of 
inflation, these primordial perturbations are sourced by scalar 
quantum vacuum fluctuations of the spacetime metric and 
inflaton field, which are first stretched to superhorizon scales 
during inflation and then reenter the horizon in the form of 
classical density perturbations after inflation. In addition to 
scalar perturbations, inflation also leads to the amplification of 
tensor perturbations of the metric, which reenter the horizon in 
the form of stochastic GWs after inflation. These primordial or 
inflationary GWs (IGWs; Grishchuk 1974; Starobinsky 1979; 
Rubakov et al. 1982; Fabbri & Pollock 1983; Abbott & 
Wise 1984) represent a prime GW signal from the early 
universe. For earlier work on the IGW interpretation of the 
signal in recent PTA data sets, see Vagnozzi (2021), 
Kuroyanagi et al. (2021), and Benetti et al. (2022). 

IGWs leave an imprint in the temperature and polarization 
anisotropies of the cosmic microwave background (CMB) 
whose relative strength compared to the contributions from 
scalar perturbations is quantified in terms of the tensor-to-scalar 
ratio, r. For the simplest type of inflation—standard single-field 
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scales, with the tensor spectral index n; being given by the so- 
called consistency relation, п, = — r/8 <0. Meanwhile, r is 
bounded from above by current CMB observations, r < 0.036 
at 9596 C.L. (Ade et al. 2021). A vanishing tensor spectral 
index, п, = 0, would imply an upper bound on the GW energy 
density spectrum of Осуй2 ~ 10 !6 at РТА frequencies, 
rendering any detection of an IGW signal in PTA observations 
hopeless. This conclusion, however, only applies to the 
standard case of single-field slow-roll inflation. Nonminimal 
scenarios may have significantly better detection prospects. 

We remain agnostic about the microphysics of inflation and 
restrict ourselves to a model-independent analysis, in which we 
parameterize the IGW signal in terms of four parameters: the 
tensor-to-scalar ratio г and tensor spectral index п, at the CMB 
pivot scale, fcwg = 0.05 Мре '/(2тао) = 7.73 х 10 1 Hz, which 
quantify the efficiency and scale dependence of GW production 
during inflation, and the reheating temperature Ти, and the number 
of e-folds during reheating №, which describe the reheating 
process after inflation. Here the factor ag in the expression for /смв 
denotes the present value of the cosmological scale factor a(t) in 
the Robertson-Walker metric; in our convention, ao — 1. 

We do not impose the standard consistency relation between r 
and n; instead, we allow both parameters to vary independently 
across large prior ranges, log, Є [—40, 0] and n; € [0, 6]. We 
note that blue values of the tensor spectral index can be generated, 
e.g., from axion-vector dynamics during inflation (Апрег & 
Sorbo 2012; Cook & Sorbo 2012; Namba et al. 2016; 
Dimastrogiovanni et al. 2017; Caldwell & Devulder 2018) or in 
other nonminimal inflation models (see Piao & Zhang 2004; 
Satoh & Soda 2008; Kobayashi et al. 2010; Endlich et al. 2013; 
Fujita et al. 2019 for an incomplete list). 

Similarly, we allow for more flexibility for Tin and М, than 
in the standard treatment of single-field slow-roll inflation. To 
illustrate this point, note that the number of e-folds during 
reheating, N,,, can be written as (Liddle & Lyth 2009) 


(14) 


1 3H24Mg, 
Na = п] 5 dI 
ЗИ + Wa) «2/3099 Т, 


where Heng is the Hubble rate at the end of inflation, ињ is the 


slow-roll inflation—the A^Qw spectrum is red-tilted at CMB equation-of-state parameter during reheating, Мы = 2.44 x 
2 
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Figure 2. Bayes factors for the model comparisons between the new-physics interpretations of the signal considered іп this work and the interpretation in terms of 
SMBHBs alone. Blue points are for the new physics alone, and red points are for the new physics in combination with the SMBHB signal. We also plot the error bars 
of all Bayes factors, which we obtain following the bootstrapping method outlined in Section 3.2. In most cases, however, these error bars are small and not visible. 
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Figure 3. Median GWB spectra produced by a subset of the new-physics models, which we construct by mapping our model parameter posterior distributions to 
I? Naw distributions at every frequency f (see Appendix B for more details and Figures 19 and 20 for the models not included here). We also show the periodogram for 
an HD-correlated free spectral process (gray violins) and the GWB spectrum produced by an astrophysical population of inspiraling SMBHBs with the parameters 
Авнв and “pup fixed at the central values ивнв of the 2D Gaussian prior distribution specified in Equation (АТ) (black dashed line). 


10/5 GeV is the reduced Planck mass, and p^ is defined below. 
We assume for definiteness that reheating is dominated by the 
coherent oscillations of the inflaton field, such that the equation of 
state is equivalent to the one of pressureless dust (i.e., matter), 
иль = 0. In typical models of single-field slow-roll inflation, one 
can often approximate Heng by the Hubble rate at the time of CMB 
horizon exit during inflation, such that 


2 2 
тае ду liu M , (15) 
3 (12/30 2 Th 
where Haaive is fixed by the tensor-to-scalar ratio г and the 
amplitude of the primordial scalar power spectrum, 
А, 2.10 x 107° (Aghanim et al. 2020), 
2 1/2 
Hena © Hnaive = (Za) Мы. (16) 


However, we already assume nonminimal dynamics in order to 
motivate a strongly blue-tilted ^ Qo spectrum, so there is no 
reason why we should make use of this approximation. In our 
analysis, we therefore treat Heng and correspondingly N, as 
independent parameters and do not fix them in terms of r and 
Та as in Equations (15) and (16). This flexibility provides us 
with more parametric freedom that we can use in order to 
ensure that the IGW signal does not violate constraints on the 
amplitude of the stochastic GWB set by the LIGO-Virgo- 
KAGRA (LVK) Collaboration (Abbott et al. 20212) and on the 
amount of dark radiation, i.e., the effective number of neutrino 
species, Мең, inferred from big bang nucleosynthesis (ВВМ) 
and the CMB (Pisanti et al. 2021; Yeh et al. 2021). As for the 


latter constraint, we specifically work with АМ. к ppn/ py, 
where ppg is the energy density of dark radiation (1.е., the 
integrated GW energy density in the context of the IGW model) 
and p, denotes the energy density of a single neutrino species. 
AN. characterizes the excess energy in radiation beyond the 
SM expectation (i.e., dark radiation) after neutrino decoupling 
and ете” annihilation, Мар = NSM + ААЛ, where ММ ~ 
3.0440 (Bennett et al. 2021). 

Under the assumptions outlined above, we are able to model 
the IGW spectrum at РТА frequencies as 


24\ g? 

Here О,/82 > 2.72 x 1075 is the current radiation energy 
density per relativistic degree of freedom, in units of the critical 
(closure) density, 80 , œ 3.93 counts the effective number of 
relativistic degrees of freedom contributing to the radiation 
entropy today, and g,(f) and gx (f) denote the effective 
numbers of relativistic degrees of freedom in the early universe 
when GWs with comoving wavenumber К = Зла reentered 
the Hubble horizon after inflation. In order to evaluate g..(f) 
апа gx (f), we use the numbers of relativistic degrees of 
freedom as functions of temperature tabulated in Saikawa & 
Shirai (2020), g,(T) and gy (Т), in combination with the 
standard temperature-frequency relation in ACDM (і.е., the 
cosmological Lambda cold dark matter SM) that follows from 
the condition К = a(T)H(T) at the time of horizon reentry. In the 
remainder of this paper, whenever we need g or gx, in a 


4/3 
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different part of our analysis, we will use the same functions 
8x(f), 8x), 8» (D), and ё„ Т). 

The inflationary dynamics give rise to the primordial tensor 
power spectrum 7, while the transfer function 7 accounts for 
the redshifting behavior of GWs after horizon reentry. In our 
analysis, we assume a constant tensor spectral tilt (1.е., zero 
running of n;) from CMB to PTA frequencies, such that 


Ff) =r As (18) 


CMB 


Meanwhile, the only relevant contribution to 7 in the PTA 
band corresponds to the transfer function that connects the 
radiation-dominated era to reheating, 


Tf) " O( fena = f 


ы OVa D І 
1 0.22071) + 0.65Cf fu. ? С 


Here the fit function in ће denominator of this expression is 
taken from Kuroyanagi et al. (2015, 2021) and describes the 
spectral turnover, Ё" — f"-?, at frequencies around f~ fin, 


f [Hz] 


Figure 4. Same as Figure 3, but for a different selection of models and showing a larger frequency range. The solid lines represent median GWB spectra for a subset of 
new-physics models (see Appendix B for more details), the gray violins correspond to the posteriors of an HD-correlated free spectral reconstruction of the 
NANOGrav signal, and the shaded regions indicate the power-law-integrated sensitivity (Thrane & Romano 2013) of various existing and planned GW interferometer 
experiments: LISA (Amaro-Seoane et al. 2017), DECIGO (Kawamura et al. 2011), BBO (Crowder & Cornish 2005), Einstein Telescope (ET; Punturo et al. 2010), 
Cosmic Explorer (CE; Reitze et al. 2019), the HLVK detector network (consisting of aLIGO in Hanford and Livingston (Aasi et al. 2015), aVirgo (Acernese 
et al. 2015), and KAGRA (Akutsu et al. 2019)) at design sensitivity, and the HL'V detector network during the third observing run (O3). All sensitivity curves are 
normalized to a signal-to-noise ratio of unity and, for planned experiments, an observing time of 1 yr. For the HLV detector network, we use the O3 observing time. 
Different signal-to-noise thresholds pmr and observing times tops can be easily implemented by rescaling the sensitivity curves by a factor of р, / /tos,. More details on 
the construction of the sensitivity curves can be found in Schmitz (2021). We emphasize that models whose median GWB spectrum exceeds the sensitivity of existing 
experiments are not automatically ruled out. This applies, e.g., to cosmic superstrings (SUPER) and the O3 sensitivity of the HLV detector network. Typically, no 
single GWB spectrum in a given model will coincide with the median GWB spectrum, which is constructed from distributions of /^Qgw values at any given 
frequency. Therefore, if the median GWB spectrum is in conflict with existing bounds, typically only some regions in the model parameter space will be ruled out, 
while others remain viable (see, e.g., Figure 11 for the SUPER model). Finally, note that any primordial GWB signal is subject to the upper limit on the amount of dark 
radiation in Equation (23), which requires the total integrated GW energy density to remain smaller than O(10- O79) (see Section 5.1). 
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which marks the end of the reheating period, 
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Ж = (20) 


with g^ = g,(Tm) and E = (Zn) and the present-day 
CMB temperature To 2.73 К (Fixsen 2009). The Heaviside 
theta function in Equation (19) denotes the endpoint of the 
IGW spectrum at f= fena, which marks the end of inflation and 
hence the onset of reheating, 


о M73 
1 Sys 


2th M/3 44/3 pr1/3 
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For fixed T,,, the frequency feng is solely controlled by Heng, 
which follows from our choice of М according to 
Equation (14) (recall that we set и = 0). ш our MCMC 
analysis, we do not sample over feng, since its precise value 
does not affect the shape of the IGW spectrum in the PTA band 
and hence the quality of our fit. Instead, we constrain its 
maximally allowed value (i.e., N,,) after identifying the viable 
region in the r-n;-7,, parameter space, in order to make sure 
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that the IGW spectrum does not violate the Ne and ГУК 
bounds. The frequency fig, on the other hand, can easily fall 
into the PTA band: from Equation. (20) we have 
Хь ~ 30 nHz(Ti/1 GeV). Therefore, we sample the reheating 
temperature in a symmetric interval around Ти = 1 GeV 
extending down to temperatures relevant for BBN, 
Tggw^-1 MeV. That is, we work with the log-uniform 
prior log, (75/1 GeV) € [—3, +3]. 


5.1.2. Results and Discussion 


We now discuss the outcome of our Bayesian fit analysis. 
First, we note that the IGW model fits the NGI5 data slightly 
better than the baseline SMBHB model. This is evident from 
the Bayes factor that we find for the IGW versus SMBHB 
model comparison, В = 8.8 + 0.3 (mean value and one 
standard deviation), and simply follows from the larger 
parametric freedom of the IGW model. Both the IGW and 
SMBHB models basically correspond to a power-law approx- 
imation of the GW spectrum. However, in the case of the IGW 
model, the parameters controlling this power law are drawn 
from prior distributions that allow for a larger amplitude and 
a steeper slope of the spectrum, which improves the quality 
of the fit. Meanwhile, the combined GW spectrum from 
inflation with an additional SMBHB signal on top compared 
to the SMBHB model alone results in a Bayes factor of 
В = 7.9 + 0.3. We thus observe a slight decrease in the 
Bayes factor, which accounts for the fact that adding the 
SMBHB signal on top of the IGW signal does not improve 
the quality of fit but merely increases the prior volume 
compared to the IGW model. 

The reconstructed posterior distributions for the parameters 
of the IGW model and its IGW--SMBHB extension are shown in 
Figure 5.5! For both models, we find а strong covariance 
between the spectral index п, and the tensor-to-scalar ratio г, 
which is approximately fit by 


n; = —0.14logyor + 0.58, 


(22) 


and which can be explained as follows: the IGW interpretation 
of the PTA signal requires the primordial tensor power 
spectrum P to take values of O(10~*) at nanohertz freque- 
ncies. This requirement fixes the parameter combination 
T (fora /fomp ) in Equation (18) and thus allows us to estimate 
the coefficients in Equation (22) as 1/logio( fora //смв) = 0.14 
and 102105 /А;) Лов тА //смв) = 0.58, respectively, where 
we used fpra = 1 nHz and P, = 0.3 x 1074, 

In addition to the strong covariance between n, and r, we note 
that the posterior probabilities of both parameters exhibit a bimodal 
distribution for both IGW and IGW--SMBHB. In the 2D distributions 
of the parameter pairs (Tn, п) and (Tu, г), this bimodality is 
accompanied by an approximate reflection symmetry with respect 
to the points (logo ЛЬ/СеУ, п) ~ (—0.5, 2.75) and 
(logo Tin/GeV, logor) ~ (—0.5, — 15), respectively. These fea- 
tures of the corner plot in Figure 5 indicate that the IGW model can 
operate in two regimes: for Т > 1 GeV, the reference frequency 
Ль is larger than the frequencies in the PTA band, and the GW 
spectrum seen by NANOGrav is composed of tensor modes that 
reentered the horizon during the radiation-dominated era. For 


5! The noise in the 95% credible interval for the nrdog,yr posterior 
distribution of the IGW--SMBHB model is due to the presence of an extended 
plateau region in the posterior distribution, which renders the kernel density 
reconstruction sensitive to Poisson fluctuations in the binned MCMC data. 
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Figure 5. Reconstructed posterior distributions for the parameters of the IGW 
(blue) and IGW--SMBHB (red) models. On the diagonal of the corner plot, we 
report the 1D marginalized distributions together with the 68% Bayesian 
credible intervals (vertical lines), while the off-diagonal panels show the 68% 
(darker) and 9596 (lighter) Bayesian credible regions in the 2D posterior 
distributions. We construct all credible intervals and regions by integrating over 
the regions of highest posterior density. The gray shaded regions are strongly 
disfavored by the NG15 data, as they result in a К ratio of less than 1/10 (see 
Equation (10)). The black shaded region results in a violation of the №, bound 
in Equation (23) (see Appendix С.І and Figure 22), assuming М, = 0 (solid 
line, № = 5 (dashed line), and М, = 10 (dotted line). Figure 21 in 
Appendix С.1 shows an extended version of this plot that includes the 
SMBHB parameters Авнв and Увнв: 


Tn < 1 GeV, on the other hand, fin can be pushed below РТА 
frequencies, and the GW spectrum in the PTA band is composed 
of tensor modes that reentered the horizon during reheating after 
inflation. In the first case, the tilt of the spectrum is directly given 
by п; in the second case, it corresponds to п, — 2. Clearly, the 
mirror symmetry in the 2D distributions of (Tin, п,) and (Та, r) is 
not exact. At the level of the GW spectrum, it is explicitly broken 
by the frequency dependence of g, and g,, as well as by the 
nontrivial shape of the transfer function in Equation (19). 

At small or large values of the reheating temperature, the 
posterior distributions develop approximately flat directions 
along the Tın axis at (n;, logiyr) ~ (2, —10) in the large-T;, 
regime and at (п,, logor) ~ (4, —20) in the small-7,, regime. 
This behavior is broadly consistent with the linear fit in 
Equation (22), the reflection symmetry discussed above, and 
the fact that, past a certain point, raising or lowering Т no 
longer influences the shape of the GW signal in the PTA band. 
A tensor index of п, = 2 at large Tin, moreover, corresponds to 
an index у = 3 in the timing residual PSD, which is among the 
best-fitting values—see the 2D posterior for A and y in 
Figure 1. The same conclusion holds at low Ти, where у = 3 is 
mapped onto n, — 2 — 2. 

Finally, we derive the Ме and ГУК bounds on the IGW 
parameter space. For Ше Мк bound, the GW spectrum, 
integrated from BBN scales to the cutoff frequency fena, must 
not exceed a certain upper limit that is set by the allowed 
amount of extra relativistic degrees of freedom at the time of 
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BBN and recombination, 


faa df 
ie 


BBN 


тов (Р) S 5.6 x 10-6 AN, (23) 
Here ввм ~ 10 '* Hz refers to tensor modes that reentered the 
Hubble horizon around the onset of BBN at T~ 10 ^ GeV 
(Caprini & Figueroa 2018), and АМ” ~ few х 0.1 denotes 
the upper limit on the amount of dark radiation. For 
definiteness, we set /ввм = 107! Hz and А HR. = 0.5 in our 
analysis; the precise АМ” value at 95% confidence level 
varies across different combinations of data sets (Pisanti et al. 
2021; Yeh et al. 2021). For given values of the parameters Tj. 
г, and п, we can then ask whether there is a cutoff frequency 
fa that leads to the saturation of the inequality іп 
Equation (23). If this is the case, f£." yields an upper bound 
on the allowed number of e-folds during reheating, Мех, 
according to Equations (14) and (21). A second constraint on 
Nn follows from the ГУК bound on the amplitude of the 


GWB (Abbott et al. 20212), 


Осу € L7 x 108 at f, = 25 Hz, Q4) 


assuming a flat GW spectrum. For a power-law spectrum, 
Осу = Осу, ЛЛук with а =n, — 2 at f 2» fm, this bound 
can be generalized following (Kuroyanagi et al. 2015, 2021) 


= 2a \!/2 Es 
Осу E TAX (> 22) Е =) , (25) 
5 Лук 


which approximately holds if a < 5/2. For given values of Ти, 
r, and n, we can then evaluate Q3, at Дух and check whether it 
exceeds the ГУК bound. If so, we place an upper bound on М, 
by demanding that иг =20 Hz (the lower end of ће LVK 
frequency band). 

The outcome of our analysis is shown Figure 22 in 
Appendix C.1. In both plots of this figure, the dependence on 
the tensor-to-scalar ratio r is eliminated by means of the linear 
relation in Equation (22). In the left panel of Figure 22, we 
present contour lines ої Му. derived from the №, bound and 
the LVK bound, respectively. In a realization of the IGW model 
with a given duration of reheating, these contour lines can be 
thought of as bounds on the 7,4—n, parameter plane—for fixed 
Nn, the regions with Му. < Nn are excluded. We find that the 
Ме bound rules out most values of the spectral index п, in the 
large-7,4 regime. At the same time, large regions of parameter 
space remain viable as long as М. < №. In fact, away from 
the region where Np” turns negative, our upper bound is 
typically rather large, № ~ O(10 --- 100), and hence easy 
to satisfy in realistic models of reheating, where 
Nn ~ ОСТ... 10). In the right panel of Figure 22, we compare 
our result to the naive expectation Ме? in single-field slow- 
roll inflation with a nearly constant Hubble rate (see 
Equation (15)). Across large parts of parameter space, we find 
that МЕ" assumes unrealistically large values, М >> 10. 

In Figure 5, we highlight the constraints on Tiy and п, (as 
well as T, and r) that we deduce from the Ме bound assuming 
Na values of М, = 0, 5, and 10. Here the constraints for 
№ = О correspond to the assumption of instantaneous 
reheating after inflation and hence represent the most 
conservative bound on the Тһ-п, parameter plane. A longer 
duration of reheating results in tighter constraints on Ти and n, 
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as illustrated by the contours for М = 5 and 10. For an even 
larger number of e-folds during reheating, see Figure 22 in 
Appendix C.1. 

In view of Figures 5 and 22, we conclude that the IGW model 
is indeed capable of fitting the NANOGrav signal across large 
regions in parameter space. An interesting viable scenario 
consists, e.g., of a large tensor spectral index, п, ~ 3 .-. 4, a tiny 
tensor-to-scalar ratio, rw 1072319, a low reheating temper- 
ature, Tj ~ 1079779) GeV, and a moderate number of e-folds 
during reheating, № < 20. It remains to be seen whether it is 
possible to identify explicit microscopic models that realize 
inflation in this parametric regime. 


5.2. Scalar-induced Gravitational Waves 
5.2.1. Model Description 


The amplitude of the primordial scalar power spectrum is 
well measured by CMB observations, A, ~ 2.10 x 10^? at the 
CMB pivot scale Ксмв = 0.05 Mpc! (Aghanim et al. 2020). If 
we naively extrapolate this value down to smaller scales, 
assuming a fixed and slightly red-tilted №Осу spectrum with 
index n,~ 0.96, we are led to conclude that there must be 
increasingly less power in scalar perturbations on shorter 
scales. This conclusion can, however, be easily avoided in 
models that deviate from the standard picture of single-field 
slow-roll inflation giving rise to a nearly scale-invariant 
spectrum of scalar perturbations. A prominent example, among 
many other mechanisms, consists in a stage of inflation close to 
an inflection point in the scalar potential, which readily 
amplifies the scalar perturbations leaving the horizon (see, 
e.g., Garcia-Bellido & Ruiz Morales 2017; Ballesteros & 
Taoso 2018; Ezquiaga et al. 2018). An enhanced scalar power 
spectrum at small scales is, therefore, a viable possibility. 
Moreover, it promises a rich phenomenology with regard to the 
production of GWs and potentially the origin of primordial 
black holes (PBHs; Carr et al. 2016; Garcia-Bellido et al. 2016; 
Inomata et al. 2017a; Inomata & Nakama 2019; Wang et al. 
2019; Escrivà et al. 2022b). The possibility of having PBH 
formation in models of single-field inflation is the subject of 
ongoing debate (Kristiano & Yokoyama 2022, 2023; Choudh- 
ury et al. 2023a, 2023b, 2023c; Firouzjahi & Riotto 2023; 
Riotto 2023a, 2023b). Below, we comment on the implications 
of this debate for our PBH-related parameter bounds. 

In cosmological perturbation theory, scalar and tensor 
perturbations evolve independently at linear order. Starting at 
second order, however, they are coupled, which means that 
large first-order scalar perturbations can act as a source of 
second-order tensor perturbations. We refer to these tensor 
perturbations, which are produced as soon as the enhanced 
scalar perturbations reenter the Hubble horizon after inflation, 
as SIGWs (Matarrese et al. 1993, 1994, 1998; Mollerach et al. 
2004; Ananda et al. 2007; Baumann et al. 2007; see 
Doménech 2021 for a review and more details on the history 
of SIGWs). At the same time, large overdensities in the tail of 
the distribution of scalar perturbations can collapse into PBHs 
upon horizon reentry. This PBH production mechanism thus 
results in a PBH population whose properties are strongly 
correlated with the spectral shape of the SIGW signal, as both 
phenomena are controlled by the scalar spectrum generated 
during inflation (Yuan & Huang 20212). For earlier works on 
the PBH/SIGW interpretation of the signal in recent PTA data 
sets, see Vaskonen & Уеегтае (2021), De Luca et al. (2021), 
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and Kohri & Terada (2021). For earlier Bayesian searches for SIGW-BOX: Flat and continuous spectral feature in Pg modeled 
an SIGW signal in PTA data sets, see Chen et al. (2020), Bian by a box function in logarithmic К space, 
et al. (2021), Zhao & Wang (2022), Yi & Fei (2023), and 
Dandoy et al. (2023), and for a search in LVK data, see Prk) = A O(ln kmax — In k)O(Ink — In К). (31) 
Romero-Rodriguez et al. (2022). 

In our analysis, we consider SIGWs in the PTA band and 
model the associated GW spectrum as follows: 


Note that the Gaussian power spectrum in logarithmic k 
space that we assume in the SIGW-GAUSS model corresponds to 
a lognormal power spectrum in linear k space. As evident from 


о \4/3 the above expressions, SIGW-DELTA represents а two-para- 

qind =Q gU) 8%. Qin 26 meter model, while SIGW-GAUSS and SIGW-BOX are three- 
ew(f) = 9,| = см C), (26) | | | 

8, B ОЛ) parameter models. Our prior choices for the respective 

parameters are listed in Table 3 in Appendix A, where we 

where the first three factors account for the cosmological use again f= k/ (27а0) to convert from wavenumber to 

redshift as in Equation (17) and the last factor denotes the GW frequency. For a given set of parameter values, we are then able 


to use the scalar power spectrum in Equation (29), 
Equation (30), or Equation (31) to evaluate the integrals in 
Equation (27) and compute the GW spectrum. For SIGW- 
= ind oo l+v GAUSS and SIGW-BOX, the integration needs to be carried out 
Qew(f) =Í dv Г du Ки, v)Pr(uk)Pr(vk). (27) numerically; for SIGW-DELTA, we can resort to the exact 
Рем analytical expression provided in Wang et al. (2019) and Yuan 

The present-day GW frequency is related to the comoving & Huang (20214). 
wavenumber by f= К/(27ао), Pr denotes the primordial 
spectrum of the comoving curvature perturbation 7, and the 
integration kernel K is given by (Espinosa et al. 2018; Кови & We now turn to the outcome of our Bayesian fit analysis. 
Terada 2018; Pi & Sasaki 2020; Gong 2022) Compared to the IGW model discussed in Section 5.1, we 


spectrum at the time of production, which we assume to be 
during the radiation-dominated era, 


5.2.2. Results and Discussion 


|.3(A2 — (1 + v? — и2)2)2 (и? + у? — 3)4 


Ka, v 
Бы 1024 uy 
2 
3 – (и + у)? 4uv " 
x |] In + лт?Ө(и + v — 43) |. (28) 
| 3 - (и — у)? и? ру? - 3 
The expression in Equation (27) illustrates the dependence of obtain even larger Bayes factors, indicating that SIGWs tend to 
the SIGW signal on the scalar input spectrum; in particular, it provide an even better fit of the NGI5 data than IGWs. 
shows that 01%, scales as 04, oc P2. We stress that the Specifically, we obtain B = 44 +2, B = 57 + 3, and 
expression in Equation (27) assumes Gaussian perturbations, B- 21 X 1 for SIGW-DELTA, SIGW-GAUSS, and SIGW-BOX, 
whose statistics are fully described by the power spectrum Pp. respectively, and B = 38 + 2, В = 59 + 3, and B= 23 + 1 
We do not consider any non-Gaussian contributions to the for SIGW-DELTA+SMBHB, SIGW-GAUSS+SMBHB, and SIGW- 


primordial scalar power spectrum in our analysis. The impact ВОХ+ӘМВНВ, respectively (see Figure 2). This improvement in 
of non-Gaussianities on the SIGW signal was studied in the quality of the fit reflects the fact that the SIGW spectra 


Nakama et al. (2017), Cai et al. (2019), Unal (2019), Atal et al. deviate from a pure power law and thus manage to provide a 


_ better fit across the whole frequency range probed by 
(2019), Yuan & Huang (2021b), Atal & Doménech (2021), й 
Adshead et al. (2021), and Ferrante et al. (2023). NANOGrav (see Figure 3). For SIGW-DELTA, we observe 


Е : à ч again that adding the SMBHB signal does not improve the 

To remain as model-independent as possible, we refrain from om of ~ = d of es volume-ot SI Eun DETA 
choosing a particular inflation model capable of generating an guay : ger p Me 

enhanced spectrum Pk. Instead, we ignore the microphysics of +$МВНВ compared to SIGW-DELTA therefore results in a slight 

inflation and work with three characteristic templates for Pr decrease of the Bayes factor. For the other two SIGW models, 

that reflect the range of possibilities that one may expect in the Bayes factors remain more or less the same, within the 


realistic models. Specifically, we consider the following statistical uncertaintity of our bootstrapping analysis. 
The reconstructed posterior distributions for all SIGW 


o ш Sharp spectral feature in Pe modeled by a models under consideration are shown in Figures 6 and 7. 
Dirac delta function in logarithmic k space, Our first conclusion in view of these results is that a successful 
explanation of the NANOGrav signal in terms of SIGWs 

Prk) = A (lnk — ША»). (29) requires the primordial scalar power spectrum to have a large 


amplitude A. We can quantify this statement in terms of the 
lower limits of the 95% Bayesian credible intervals for A that 
we obtain from integrating the corresponding 1D marginalized 
3 posteriors over the regions of highest posterior density: for 
ЕЕ - mts) | (30) SIGW-DELTA, SIGW-GAUSS, and SIGW-BOX, we respectively 

A find log,,A 2, 4.55, log A Z 4.43, and log, А Z 4.90. At 


SIGW-GAUSS: Broad spectral feature in Pr modeled Бу а 
Gaussian peak in logarithmic К space, 


Prk) = 


A 
От A 
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the same time, the enhanced power in scalar perturbations must 
be localized at the right scales, so that the resulting SIGW 
signal falls into the PTA band. This requirement leads to 
bounds on the frequencies fhin» fina and №» that can be read off 
from Figures 6 and 7 and that are summarized in Table 4. 

A notable feature in this context is that the posterior 
distributions for fj Лау, and f all extend to large 
frequencies, much beyond the upper limit of the NANOGrav 
band. The reason for this is that the NGI5 data are best fit by 
the low-frequency tail of the SIGW spectrum (see Figure 19). 
Beyond the NANOGrav band, the SIGW spectrum keeps 
increasing until it reaches its maximal value at f>> 1 nHz. This 
observation also explains the flat directions in the 2D posterior 
distribution for A and fain in Figure 6 and the 2D posterior 
distributions for A and f, in Figure 7. A simultaneous increase 
in A and fain orf along these flat directions moves the peak in 
the GW spectrum to higher frequencies, but it preserves the 
shape of the low-frequency tail in the PTA band and hence 
does not affect the quality of the fit. 

We also observe that the 2D posterior distributions for A and 
Лу for SIGW-DELTA and SIGW-GAUSS are roughly similar to 
each other, with the distribution for the SIGW-GAUSS model 
being slightly broader. This result is consistent with the fact 
that SIGW-DELTA and SIGW-GAUSS are nested models; SIGW- 
DELTA can be obtained from SIGW-GAUSS in the limit A — 0. 
The slightly broader posterior distribution for A and f in the 
right panel of Figure 7 thus reflects the extra dimension in the 
parameter space of the SIGW-GAUSS model and the additional 
parametric freedom that comes with it. 

Finally, we comment on the bounds on the parameter space 
of the three SIGW models. As in the case of the IGW model, we 
identify regions of the parameter space where the K ratio in 
Equation (10) falls below 1/10. In these regions, which are 
shaded in gray in Figure 7 and labeled №015, adding the SIGW 
contribution to the GW signal leads to much worse agreement 
with the data than in the case of no SIGW contribution at all. In 
fact, parameter points in these regions lead to an SIGW signal 
that exceeds the observed signal—in other words, the gray 
shaded regions are ruled out because they result in too strong of 
an SIGW signal. For SIGW-BOX, we are not able to place a K- 
ratio bound on the 2D parameter space spanned by A and fi, 
because we additionally marginalize over / ах» That is, for any 
pair of values of A and Дід» the SIGW signal can be made 
arbitrarily small if we choose f,,,, close enough to / 

These bounds are an important result of our analysis that is 
independent of the origin of the NANOGrav signal. They 
provide valid constraints on the parameter space for both the 
SIGW-DELTA and SIGW-GAUSS models, regardless of whether 
these models contribute to the explanation of the observed 
GWB. In addition, they are complementary to other existing 
bounds, in particular, the requirement that the mass density of 
PBHs produced alongside SIGWs must not exceed the energy 
density of DM, 

вн < 1. (32) 
where jenu = Орвн/Орм denotes the PBH DM fraction 
integrated over the entire PBH mass spectrum. We provide 
more details on how we compute јьвн in Appendix C.2; here 
we simply report our final results in terms of the teal contour 
lines labeled вн = 1 in Figures 6 and 7. For SIGW-BOX, the 
PBH bound in the А- f, plane is computed for fin =107!! Hz 


(the lower end of our prior range), and the PBH bound in the 
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Figure 6. Same as in Figure 5, but for the SIGW-BOX model. The regions above 
the teal contour lines labeled /рвн = 1 lead to the overproduction of PBHs, 
according to our analysis in Appendix C.2; however, see text for more 
discussion. Figure 23 in Appendix C.2 shows an extended version of this plot 
that includes the SMBHB parameters Авнв and Эвнв- 


10810 A 


А—/ и plane is computed for f,,,,=10~> Hz (the upper end of 
our prior range). For SIGW-GAUSS, we show the PBH bound in 
the A-f, plane for А = 1 (solid contour line) and A=2 
(dashed contour line). 

In all three cases, we find that the PBH bound is very 
restrictive, ruling out most of the parameter space favored by 
the NGI5 data. If taken at face value, the PBH bound therefore 
renders the SIGW interpretation of the NANOGrav signal less 
likely. However, we stress that the computation of /рвн is very 
sensitive to various assumptions and numerical steps in the 
analysis. Slight changes in the computational strategy may 
result in very different results for /рвн, which is why the results 
reported in Figures 6 and 7 need to be treated with caution. In 
view of the large conceptional uncertainty in the computation 
Of /рвн, one needs to be careful not to draw any premature 
conclusions. At the same time, the PBH bounds in Figures 6 
and 7 illustrate that small regions of parameter space do remain 
viable. In fact, for SIGW-DELTA and SIGW-GAUSS, it is even 
possible to realize />вн = 1 inside the 68% credible regions. 
That is, along the /рвн = 1 contour lines and inside the 68% 
credible regions, we find scenarios where SIGWs manage to 
provide a good fit to the NGI5 data and PBHs account for the 
entire DM in our universe. 

In closing, we remark that the above conclusions are endangered 
by the recent claim of a no-go theorem for PBH formation from 
single-field inflation (Kristiano & Yokoyama 2022, 2023). Kris- 
tiano & Yokoyama (2022, 2023) argue that enhanced scalar 
perturbations at small scales lead to unacceptably large one-loop 
corrections to the scalar power spectrum at large scales. In terms of 
the model parameters discussed in this section, this means that the 
amplitude A must be small, log, A < -2 otherwise, the loop 
corrections to the scalar power spectrum will exceed the measured 
amplitude at CMB scales, А, 2.10 x 107°. At present, this claim 
is subject to an ongoing debate; it is notably challenged by Riotto 
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Figure 7. Same as in Figure 6, but for the SIGW-DELTA (left panel) and SIGW-GAUSS (right panel) models. The solid and dashed teal contour lines labeled вн = 1 in 
the right panel indicate the PBH bound for A — 1 and A — 2, respectively. Figure 24 in Appendix C.2 shows extended versions of the two plots that include the 


SMBHB parameters Авнв and 7внв: 


(2023a, 2023b) and Firouzjahi & Riotto (2023). However, if it 
should prove to be valid, the requirement of log; A < —2 would 
clash with the lower bounds on A listed above. In this case, one 
would then have to either give up on the SIGW interpretation of 
the NANOGrav signal or seek to construct inflation models that 
evade the arguments of Kristiano & Yokoyama (2022, 2023) and 
still lead to a large SIGW signal. 


5.3. Cosmological Phase Transitions 
5.3.1. Model Description 


In the cosmological context, first-order phase transitions take 
place when a potential barrier separates the true and false 
minima of a scalar field potential." The phase transition is 
triggered by quantum or thermal fluctuations that cause the 
scalar field to tunnel through or fluctuate over the barrier, 
which results in the nucleation of bubbles within which the 
scalar field is in the true vacuum configuration. If sufficiently 
large, these bubbles expand in the surrounding plasma where 
the scalar field still resides in the false vacuum. The expansion 
and collision of these bubbles (Kosowsky et al. 1992a, 1992b; 
Kosowsky & Turner 1993; Kamionkowski et al. 1994; Caprini 
et al. 2008; Huber & Konstandin 2008), together with sound 
waves generated in the plasma (Giblin & Mertens 2013, 2014; 
Hindmarsh et al. 2014, 2015), can source a primordial GWB 
(see Witten 1984; Hogan 1986 for seminal work). For earlier 
Bayesian searches for a phase transition signal in PTA data, see 
Arzoumanian et al. (2021) and Xue et al. (2021). 


82 The scalar field can either be an elementary field of a weakly coupled theory 
or correspond to the vacuum condensate of a strongly coupled theory. 
Scenarios with several scalars are also possible. 

83 Turbulent motion of the plasma can also source GWs; however, its 
contribution is usually subleading compared to the two other contributions (see 
the discussion in Caprini et al. 2020). Therefore, we ignore GWs sourced by 
turbulence in our analysis. 
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Generally, the GWB produced during the phase transition is 
a superposition of the bubble and sound-wave contributions. 
However, if the bubble walls interact with the surrounding 
plasma, most of the energy released in the phase transition is 
expected to be converted to plasma motion, causing the sound- 
wave contribution to dominate the resulting GWB. An 
exception to this scenario is provided by models in which 
there are no (or only very feeble) interactions between the 
bubble walls and the plasma, or by models where the energy 
released in the phase transition is large enough that the friction 
exerted by the plasma is not enough to stop the walls from 
accelerating (runaway scenario). However, determining 
whether or not the runaway regime is reached is either model 
dependent or affected by large theoretical uncertainties (see, 
e.g., Ai et al. 2023; Krajewski et al. 2023; Li et al. 2023, for 
recent work on this topic). Therefore, in this work, we perform 
two separate analyses: a sound-wave-only analysis (РТ- 
SOUND), where we assume that the runaway regime is not 
reached and sound waves dominate the GW spectrum, and a 
bubble-collisions-only analysis  (PT-BUBBLE), where ме 
assume that the runaway regime is reached and bubble 
collisions dominate the GW spectrum. 

We parameterize the GWB produced by both sound waves 
and bubble collisions in a model-independent way in terms of 
the following phase transition parameters: 


1. T,, the percolation temperature, i.e., the temperature of 
the universe when ~34% of its volume has been 
converted to the true vacuum (Ellis et al. 2019a). For 
weak transitions, this temperature coincides with the 
temperature at the time of bubble nucleation, T, ~ Ту. 
Conversely, for supercooled transitions, we typically 
have T,« Tą. Barring the case of extremely strong 
transitions, ах 1 (see below), which we do not 
consider in this work, Т, also determines the reheating 
temperature after percolation, Ти ~ Т. (Ellis et al. 
20192). 


THE ASTROPHYSICAL JOURNAL LETTERS, 951:L11 (56pp), 2023 July 1 


ШИШЕ rpT-BUBBLE + SMBHB 
ШШШ PT-BUBBLE 


ж o 
Е 
a0 
$-1 
cé 
gl 
= 
2 
8 2 
nana nala ПРОЧАНИ О о алаа 


iii 
—1 0 


logio T. /GeV 10510 @* 


Afzal et al. 


EN pT-sOUND + SMBHB 
ШШШ PT-sOUND 


I 2,3, | ша 
-2 -1 
10510 Н, Е» 


10510 T. /GeV 


19510 Ox 


Figure 8. Same as in Figure 5, but for the PT-BUBBLE (left panel) and PT-SOUND model (right panel). Figure 25 in Appendix C.3 shows the same plots but with the 
parameter a fixed by causality, a — 3. Figures 26 and 27 in Appendix C.3 show extended versions of the two plots that include the spectral shape parameters a, b, c 


and the SMBHB parameters AguB and Эвнв. 


2. o, the strength of the transition, 1.е., the ratio of the 
change in the trace of the energy-momentum tensor 
across the transition and the radiation energy density at 
percolation (Ellis et al. 2019b; Caprini et al. 2020). 

3. Н.К = В/Н, 21, the average bubble separation in units 
of the Hubble radius at percolation, Ну 1 For relativistic 
bubble velocities, which is what we consider in the 
following, R, is related to the bubble nucleation rate, д, 
by the relation HR, = (87) ЗН, / B. 


In addition to the parameters Tx, ах, and H,R,, the GWB 
produced by a phase transition also depends on the velocity of 
the expanding bubble walls, v,,. However, deriving the precise 
value of this quantity is an open theoretical problem, which 
depends on model-dependent quantities, such as the strength of 
the interactions between the bubble walls and the SM plasma. 
Therefore, in our analysis, we fix the bubble velocity to unity 
(1.е., the speed of light in natural units). This assumption is well 
justified for strong phase transitions (Bodeker & Moore 2017), 
which, realistically, are the only ones that could lead to a 
detectable signal in our current data. In particular, we fix v,, 
= 1 for both phase transition scenarios that we are interested in, 
PT-SOUND and PT-BUBBLE. In the latter case, v, — 1 is 
automatically implied by the runaway behavior of the phase 
transition; in the former case, one actually expects a subluminal 
terminal velocity, v,, < 1. In this sense, our decision to fix у, 
= | amounts to the optimistic assumption that this terminal 
velocity is numerically close to v,, — 1. A similar approach is 
followed by the authors of the LISA review paper (Caprini 
et al. 2020), who work with v,, — 0.95 throughout most of their 
analysis in the absence of more detailed microphysical 
calculations. Finally, we point out that the parameterization 
of the GWB signal in terms of H,R, = (81)? v, Ну, / 0 
already absorbs a large part of the dependence on the bubble 
wall velocity. The remaining v, dependence is mostly 
contained in the efficiency factor к, (see below). However, in 
the regime of large a, values, o ~ 0.3 --- 10, which turn out 
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to be preferred by the NGI5 data (see Figure 8), this 
dependence is rather weak (see Figure 13 in Espinosa et al. 
2010), which justifies again keeping v,, fixed. 

The GWB spectrum sourced by bubbles and sound waves 
can be written in terms of these parameters as 


2 
90Р = ом = ) (IRA SCf /},) (33) 
1 + Ax 
^ Ky Q т 
Qf) =Р axe, A} RASC /f.). (34) 
1 + Ax 


Here €, = 0.0049 (Jinno & Takimoto 2017) and 0, = 0.036 
(Hindmarsh et а. 2017) the efficiency factor 
к; = аж/(0.73 + 0.083 /ағ + ах) (Espinosa et al. 2010) 
gives the fraction of the released energy that is transferred to 
plasma motion in the form of sound waves, and D accounts for 
the redshift of the GW energy density, 


eq \4/3 

2 TÉ E 

Е g —® ~ 1.67 x 107—5, (35) 
90 Mg Hj Ey, 


We recall that То and Но denote the photon temperature and 
Hubble rate today. The degrees of freedom gą and gx, in 
Equation (35) are evaluated at T = Ту, and о, is the number of 
degrees of freedom contributing to the radiation entropy at the 
time of matter-radiation equality. The production of GWs from 
sound waves stops after a period 7,,,, when the plasma motion 
turns turbulent (Ellis et al. 2019a, 2019b, 2020; Guo et al. 
2021). In Equation (34), this effect is taken into account by the 
suppression factor 


Trw) = 1 (1 + 25А) 7^, (36) 
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where the shock formation timescale, 7,,,, can be written in 
terms of the phase transition parameters as Tew ^: Rx / Ор, with 
О? яз 35,04 /[4(1 + a)] (Weir 2018). 

The functions 55, describe the spectral shape of the GWB 
and are expected to peak at the frequencies 


eq \1/3 

7 [53 d т, 
я dis 1 GeV 
where the values of the peak frequencies at the time of GW 
emission are given by f; —0.58/R, (Jinno & Takimoto 2017) 
and f*=1.58/R, (Hindmarsh et al. 2017). In passing, we 
mention that the numerical factors in these estimates may still 
change in the future, as our understanding of cosmological 
phase transitions improves. However, at the level of our 
Bayesian fit analysis, changes in these prefactors can be 
absorbed in the temperature scale Tą, which in its role as an 
independent fit parameter only controls the peak frequencies in 
Equation (37). A similar argument applies to the numerical 
factors in Equations (33) and (34): changes in these prefactors 
can always be absorbed in a rescaled version of the fit 
parameter a,, which only appears in the expressions for the 
peak amplitudes of the GWB signal. 


The shape of the spectral functions can be usually 
approximated with a broken power law of the form 


1 — (a4 by 
N (Ьх—“/© 15 ах?/сус` 


fy, Rs 
HRe 


$, = 48.5 nHz (37) 


Sx) = (38) 


Here a and b are two real and positive numbers that give the 
slope of the spectrum at low and high frequencies, respectively; 
c parameterizes the width of the peak. The normalization 
constant, Л/, ensures that the logarithmic integral of S is 
normalized to unity and is given by 


м= (2) (=) Ta/mV/m) (39) 
a b n Г(с) 


where n=(a+b)/c and Г denotes the gamma function. 
While the values of the coefficients a, b, and c can in principle 
be estimated from numerical simulations, we allow them to 
float within the prior ranges given in Table 3. These prior 
ranges were chosen to roughly reflect the typical uncertainties 
of numerical simulations and any possible model dependency 
of these coefficients (see, e.g., Hindmarsh et al. 2017, 2021; 
Cutting et al. 2018, 2021; Lewicki & Vaskonen 2020, 
a) 


5.3.2. Results and Discussion 


The reconstructed posterior distributions for the parameters 
О, Ту, and HR, of the PT-SOUND and PT-BUBBLE models 
are reported in Figure 8, both for the case where the phase 
transition is assumed to be the only source of GWs (blue 
contours) and for the scenario where we consider the 


ga Causality fixes the spectral index of the phase transition GWB signal to 
a=3 in the low-frequency limit. However, given the simple power-law 
parameterization adopted in this work, double-peak features observed in the 
results of numerical simulations (Hindmarsh et al. 2017; Hindmarsh & 
Hijazi 2019) can appear as a deviation from this behavior near the peak 
frequency. Nonetheless, in Appendix C.3, we also report the results of an 
analysis in which the low-frequency slope is fixed to a = 3. 
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superposition of the phase transition and SMBHB signals 
(red contours)? Corner plots including the posterior distribu- 
tions for the spectral shape parameters a, b, c and SMBHB 
parameters Авнв and Увнв are reported in Figures 26 and 27 in 
Appendix C.3. 

In all analyses, the data prefer a relatively strong and slow 
phase transition. Specifically, for PT-BUBBLE, we find a, > 1.1 
(0.29) and H.R, > 0.28 (0.14) at the 68% (95%) credible level. 
When the SMBHB signal is added on top of the GWB 
predicted by PT-BUBBLE, we find а,» 1.0 (0.23) and 
HR, > 0.26 (0.11) at the 68% (95%) credible level. For the 
PT-SOUND model, ме find 40,720.42 (0.37) апа 
HR, Є [0.053, 0.27] ([0.046, 0.89]) at the 68% (95%) credible 
level. Including the SMBHB signal on top of the one predicted 
by PT-SOUND, we find а є [0.46, 5.4] (20.16) and 
Н.К, € [0.054, 0.35] (20.0015) at the 68% (95%) credible 
level. 

As can be seen in Figure 3, for both phase transition models, 
the reconstructed GWB spectrum tends to peak around the 
higher frequency bins and fit the signal in the lower frequency 
bins with the left tail of the spectrum. Specifically, for the PT- 
BUBBLE model we find Tẹ Є [0.047, 0.41) ([0.023, 1.75]) GeV 
at the 6896 (9596) credible level, whereas for the PT-SOUND 
model we get Т, Є [4.7, 33] ([2.7, 93]) MeV at the 68% (95%) 
credible level. The shift between these Tą intervals is partially 
explained by the different numerical factors in the frequencies 
/Ё апа i, (see Equation (37)). As explained below 
Equation (37), any change in these numerical factors can be 
reabsorbed in a redefinition of the fit parameter Ту. 

The inclusion of the SMBHB signal, by adding power to the 
lowest frequency bins, allows the T, posterior for the PT- 
SOUND model to extend to higher values. In this case, we find 
that T, € [4.9, 50] ([0.8, 2 x 10°]) MeV at the 68% (95%) 
credible level. Here the increase in the 68% upper limit is 
reflected in the slight shift between the red and blue dashed 
vertical lines in the 1D marginalized posterior distribution for 
T, in the right panel of Figure 8. The drastic increase in the 
95% upper limit, on the other hand, is related to the fact that 
adding the SMBHB signal to the GWB results in a flat plateau 
region in the posterior distribution of the PT-SOUND model 
parameters where the NANOGrav signal is mostly explained 
by the SMBHB contribution to the GWB. The 95% credible 
regions for the PT-SOUND--SMBHB model cover much of this 
plateau, which explains their large extent and noisy appearance 
in Figure 8. For the PT-BUBBLE model, the inclusion of the 
SMBHB signal is less significant, and we find Ту Е [0.046, 
0.46] ([0.017, 3.27]) GeV at the 6896 (9596) credible level. 

The larger phase transition temperatures observed for the PT- 
BUBBLE model are a consequence of the smaller value of the 
peak frequency at the time of emission, hu but also of the 
lower prior range for the low-frequency spectral index adopted 
for the PT-BUBBLE model. Indeed, a shallower low-frequency 
tail allows spectra with a higher peak frequency to still produce 
a sizable signal in the lowest frequency bins. In Appendix C.3, 
we report the results of the analysis in which the low-frequency 
slope is set to the value predicted by causality (a — 3). In this 
case, as expected, the reconstructed phase transition 


35 The noise in the 95% credible regions of the posterior distributions of the 
PT-SOUND+SMBHB model is due to the presence of an extended plateau region 
in the posterior distribution, which renders the kernel density reconstruction 
sensitive to Poisson fluctuations in the binned MCMC data. 
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temperatures for the two phase transition models are closer to 
each other. 

The corner plots in Figure 8 also illustrate that, as expected 
from the expression for the peak frequency in Equation (37), 
there is an approximately linear correlation between 105, Tx 
and log, А.К. For a, < 1, we instead find a, ~ 1/(H.R..) for 
the PT-BUBBLE model and 02 c 1/(H,Rx) for the PT-SOUND 
model as expected from the factors in Equations (33) and (34). 

We also notice that, for both models, the posterior distribution 
for Ту is peaked at significantly larger values compared to what 
was derived in the 12.5 yr analysis (Arzoumanian et al. 2021). 
This shift results from the reconstructed ^ Qc spectrum being 
bluer than the one derived for the common process observed in 
the 12.5 yr data set. As a result, the lowest frequency bins, which 
were fit by the high-frequency tail of the phase transition 
spectrum in the 12.5 yr analysis, are now fit by the low- 
frequency tail of the spectrum. This then translates into a higher 
peak frequency and therefore a higher transition temperature. 
Incidentally, the larger reconstructed value for the transition 
temperature allows the phase transition signal to safely evade 
bounds from BBN and CMB observations (Bai & Korwar 2022; 
Deng & Bian 2023) for both the PT-BUBBLE and PT-SOUND 
models, which constrain the phase transition parameter space at 
temperatures around Т, ~ 1 MeV. 

Instead, we conclude that the reconstructed posterior 
distribution of Ту is compatible with phase transition scenarios 
that have been discussed in the literature as a possible source of 
GWs in the PTA band: (1) BSM models in which the chiral- 
symmetry-breaking phase transition in quantum chromody- 
namics (QCD) is a strong first-order phase transition (see, e.g., 
Li et al. 2021; Neronov et al. 2021), and (ii) strong first-order 
phase transitions in a dark sector composed of new BSM degrees 
of freedom (see, e.g. Nakai et al 2021; Ratzinger & 
Schwaller 2021). ш view of the NGI5 data, both of these 
options for the particle physics origin of the phase transition 
signal remain viable. A third option may consist in a strongly 
supercooled first-order electroweak phase transition (Kobakhidze 
et al. 2017). 

Finally, we report that, like the models studied in Sections 5.1 
and 5.2, the phase transition models provide a better fit of the 
№015 data than the base SMBHB model. The Bayes factors for 
PT-BUBBLE and PT-SOUND versus SMBHB are В = 18.1 + 0.6 
and B = 3.7 + 0.1, respectively, while the Bayes factors for PT- 
BUBBLE+SMBHB and PT-SOUND+SMBHB versus SMBHB are 
В = 12.6 + 0.5 and В = 6.5 + 0.3, respectively. An interest 
ing observation in view of these results is that adding the 
SMBHB contribution to the GWB signal does not help to 
improve the quality of the fit for PT-BUBBLE—in this case, we 
find again a decrease in the Bayes factor going from PT-BUBBLE 
to PT-BUBBLE--SMBHB because of the larger prior volume—but 
it does lead to a better fit for PT-SOUND. This model benefits 
from the additional SMBHB contribution because it can add 
power to the low frequency bins in the GW spectrum that the PT- 
SOUND model struggles to fit well on its own (see Figure 3). The 
reason for this, in turn, is the prior range for the spectral index at 
low frequencies, a, which can as be as low as a— 1 for PT- 
BUBBLE, but which we require to be at least а = 3 for PT-SOUND 
(see Table 3). Another consequence of this interplay between the 
phase transition and SMBHB signals is that the NANOGrav 
signal may in fact be dominated by SMBHBs. This possibility is 
realized when the PT-SOUND model parameters fall into the 
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plateau region in Figure 8 (i.e., the red 9596 credible regions in 
the right panel) and the SMBHB parameters are close to 
logo Annus ^ —(15 --- 14) and Эвнв ~ 3 --- 4 (see Figure 27). 


5.4. Cosmic Strings 
5.4.1. Model Description 


Cosmic strings are effectively 1D topological defects that 
can form in the early universe as a consequence of a 
cosmological phase transition (Kibble 1976). Rigorously 
speaking, the criterion for the formation of a cosmic-string 
network is that the underlying phase transition must entail the 
spontaneous breaking of a local or global symmetry and end on 
a vacuum manifold with a nontrivial first homotopy group. In 
practice, the most relevant scenario satisfying this criterion is 
the cosmological breaking of a U(1) symmetry. The breaking 
of global U(1) symmetries results in the formation of global- 
string networks, which happens, e.g., in axion models. The GW 
signal in this case is suppressed because global strings lose 
most of their energy by emitting light pseudo-Nambu-Gold- 
stone bosons (1.е., *axions"). In the following, we therefore 
focus on local strings from the spontaneous breaking of a local 
U(1) symmetry as occurs in many particle physics models of 
the early universe, such as grand unified theories (GUTs), 
where intermediate symmetry breaking stages can be realized 
in the form of string-producing phase transitions. Cosmic 
strings are thus well motivated by ideas about particle physics 
at very high energies, which prompts us to consider them as yet 
another possible source of exotic GWs (Vilenkin 1981a; Hogan 
& Rees 1984; Vachaspati & Vilenkin 1985; Accetta & 
Krauss 1989; Bennett & Bouchet 1991; Caldwell & 
Allen 1992). For earlier work on the cosmic strings interpreta- 
tion of the signal in recent PTA data sets, see, e.g., Blasi et al. 
(2021), Ellis & Lewicki (2021), and Blanco-Pillado et al. 
(2021). 

In our analysis, we study ordinary stable cosmic strings, as 
well as metastable strings, which are unstable against the 
nucleation of GUT monopoles. In passing, we also comment on 
cosmic superstrings, which do not have a particle physics 
origin but are present in some string-theoretic models of the 
early universe. Our starting point for describing all these 
scenarios is the Nambu-Goto action, which treats cosmic 
strings as featureless 1D objects that can be characterized by a 
single parameter: their tension, i.e., energy per unit length, u. In 
the case of metastable strings, there is a second parameter, к, 
which controls the lifetime of the network. For superstrings, the 
relevant parameters are u and the intercommutation probability 
P (see below). Before we move on, we note that, as an 
alternative approach to the Nambu-Goto approximation, it is 
possible to describe cosmic strings in a field-theoretical 
language, e.g., in terms of the Abelian-Higgs model. The 
comparison of these two approaches, i.e., the relation between 
Nambu-Goto strings and Abelian-Higgs strings, is the subject 
of ongoing research (Blanco-Pillado et al. 2023), and we refer 
the reader to Section 3.4 of Auclair et al. (2020) for an extended 
discussion. Furthermore, we only consider the GWB produced 
by cosmic-string loops in our analysis and disregard the 
subdominant contribution from long (infinite) strings (see a 
discussion of this point in Section 4.4 of Auclair et al. 2020). 

Shortly after their formation, cosmic strings enter a scaling 
regime where the total energy stored in the network, pes, 
remains a constant fraction of the critical energy density, 
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Qes = р../р. ғ const (Hindmarsh & Kibble 1995). This 
behavior is possible because long strings, stretching over 
superhorizon distances, frequently intercommute with each 
other, thereby producing an abundance of closed string loops 
on subhorizon scales that radiate energy in the form of GWs. 
These GWs are sourced by the oscillations of the loops under 
their own tension, as well as by localized features (“cusps” and 
"Kinks") propagating along the loops. The superposition of the 
GWs emitted by the individual loops in the network thus results 
in a stochastic GWB, 


K mas 
ФР = 27 (Сп)? ЎА Py TAG. (40) 
ЗНо = 


Here the dimensionless factor Gu denotes the cosmic-string 
tension in units of Newton’s constant, for which we choose a log- 
uniform prior in our numerical analysis, logo Gu € [—14, —6] 
for stable strings and superstrings and logi; Си Є [—14, —1.5] 
for metastable strings. The sum in Equation (40) runs over the 
harmonic excitations of the closed string loops that, given a loop 
of length £, correspond to oscillations at frequency f= 2k/£. We 
evaluate the sum starting at the fundamental oscillation mode, 
k= 1, and including terms up to kmax = 10°, which ensures good 
convergence of the GW spectrum. 

The dimensionless factor P, inside the sum describes the 
GW power, in units of Gy, that is emitted by a loop in its kth 
excitation. For large kmax, we can write 


Г 1] 


=, 41 
t си i 
where the prefactor Г/С (q) ensures that the total power emitted 
in GWs amounts to >Р; = Г and where the power-law index q 
depends on the predominant source of GWs on the loops. In 
our analysis, we set Г = 50, as suggested by numerical 
simulations (Blanco-Pillado & Olum 2017), and discuss 
different possibilities for g. The actual average GW spectrum 
from non-self-intersecting loops is still uncertain. We therefore 
choose several different models that give an idea of the range of 
possibilities. Specifically, we consider four different models of 
stable cosmic strings: 

STABLE-C: Stable strings emitting GWs only in the form of 
GW bursts from cusps on closed loops, q — 4/3 (Vachaspati & 
Vilenkin 1985). 

STABLE-K: Stable strings emitting GWs only in the form of 
GW bursts from kinks on closed loops, q = 5/3 (Damour & 
Vilenkin 2001). 

STABLE-M: Stable strings emitting monochromatic GWs at 
the fundamental frequency f= 2/0 of closed loops. 

STABLE-N: Stable strings described by numerical simulations 
including GWs from cusps and kinks (Blanco-Pillado et al. 
2011, 2015). 

For STABLE-N, P, is dominated by cusp emission at large k, 
i.e., q724/3 for > 100, while at smaller k it deviates from the 
pure cusp spectrum, reaching q values of up to q ~ 1.5 around 
k~ 10. Meanwhile, Equation (41) is irrelevant for STABLE-M; 
for this model, we simply set P, -Г if k=1 and P, = 0 
otherwise. More details on our four stable-string models can be 
found in Blanco-Pillado et al. 2021. 

Finally, the frequency-dependent factor 7, in Equation (40) 
represents an integral of the number density of closed string 
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loops, пЁ, t), over all possible GW emission times, 
Р 5 
f Чы a(to) f a(to) 


Here ig is the present time and fini is the time when the network 
reaches the scaling attractor solution. The precise value of tini 
only affects the high-frequency part of the GW spectrum and 
plays no role in our analysis. The loop number density n; can be 
estimated based on the velocity-dependent one-scale (VOS) 
model, 

"ERE: Cx O(t — te) O(te — iad (2s |, (43) 
oo + ГСи + ды) La) 


where F = 0.1 is an efficiency factor (Blanco-Pillado et al. 
2014; Auclair et al. 2020) and an asterisk indicates that the 
corresponding quantity needs to be evaluated at the time of 
loop formation, which follows from solving the following 
relation for tą: 


| (44) 
ax + ГСи 


The time-dependent functions C and a characterize the 
efficiency of loop formation from the network and the typical 
loop size at the time of production, respectively, 


_ 6 70 
42 £80) 


Here & is the loop-chopping parameter and o; controls the loop 
length at the time of production in units of the correlation 
length of the long-string network, Г, = &1. We set б = 0.23 and 
ат, 0.37, in agreement with numerical simulations (Martins & 
Shellard 2002; Blanco-Pillado et al. 2011, 2014). With € fixed, 
we are able to solve the VOS equations for the dimensionless 
correlation length £ and the rms velocity of the long strings, v, 
and hence determine the time dependence of C and a. In doing 
so, we account for the exact evolution of the scale factor, a, in 
ACDM, including the temperature-dependent variation in the 
number of relativistic degrees of freedom. At very high 
temperatures, this analysis returns 6, = 0.27 and 7, = 0.66, 
such that а œ 0.10 deep in the radiation-dominated era. The 
loop number density n; obtained in this way agrees very well 
with the result of numerical simulations (Blanco-Pillado et al. 
2014) in the limit of constant g, and gx s- 

In the case of metastable strings, the loop number density in 
Equation (43) receives two correction factors, 


C(t) a(t) = ar €t). (45) 


п" (р, t) = Oft, - tx) E(L, Dni, t). (46) 


The Heaviside function accounts for the fact that we expect 
loop formation to cease when monopole nucleation becomes 
efficient and the network transitions from the scaling regime to 
the decay regime at times around 


1 


= p (47) 
" 2 
Ги 
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with the decay rate, Го, counting the number of monopole 
nucleation events per time and per unit string length, 


(48) 


Неге к is a measure for the ratio of the GUT and U(1) 
symmetry breaking scales in the underlying GUT model. 
Specifically, /K describes the ratio of the GUT monopole mass 
and the square root of the U(1) string tension, 


MGUT 
1/2 
H / 


AGUT 
- GUT 


Vi = 
Ava) 


(49) 


Meanwhile, the second new factor in Equation (46) represents 
an exponential suppression term, reflecting the depletion of the 
loop number density in the decaying network, 


‚1 2 
EM, се що гона? |. 


As evident from Equations (48) and (50), the loop number 
density of a metastable network is exponentially sensitive to the 
decay parameter к. For Vk >> 10, the lifetime of the network 
exceeds the age of the universe, such that the resulting GW 
signal is indistinguishable from the signal of a stable-string 
network. For va ~ 1, on the other hand, the network decays 
very fast in the early universe, such that no GW signal at PTA 
frequencies is generated. For these reasons, we choose a 
uniform prior on /& in a rather narrow range, Vk Є [7.0, 9.5], 
which is enough to capture the nontrivial aspects of the 
metastable scenario. 

Monopole nucleation in a metastable network results in 
string segments, dumbbell-shaped pieces of string with 
monopoles attached on either end. In many GUT models, 
these monopoles still carry unconfined flux, such that we 
actually need to distinguish between monopoles and antimono- 
poles. Monopoles and antimonopoles with unconfined flux 
occur, e.g., in GUT models involving, schematically, the 
symmetry breaking pattern 


SU (2) x СА» ^ UC), x U(1)5 ^ UWI). 


(50) 


(51) 


In this case, we expect the string segments to dissipate most of 
their energy via the emission of massless vector bosons soon 
after their formation. Alternatively, monopoles may carry no 
unconfined flux, in which case they are able to lose energy only 
via the emission of GWs. In our analysis, we cover both 
possibilities: 

META-L: Metastable strings; monopoles carry unconfined 
flux; GW emission from loops only. 

META-LS: Metastable strings; monopoles carry no uncon- 
fined flux; GW emission from loops and segments. 

For META-LS, we also require the number density of 
segments that result from the decaying network; the relevant 
expressions are summarized in Appendix C.4. Moreover, we 
need to specify the power spectrum of GWs emitted by string 
segments. Following Buchmuller et al. (2021), we use again 
Equation (41) and set q — 1, which can be analytically derived 
in the straight-string approximation (Martin & Vilenkin 1997). 
At the same time, we assume a pure cusp spectrum (q — 4/3) 
for the GW emission from loops in both META-L and META-LS. 

Finally, we also consider cosmic superstrings: 

SUPER: Cosmic superstrings; suppressed intercommutation 
probability P; GWs from cusps on loops, q = 4/3. 
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Cosmic superstrings are characterized by a reduced inter- 
commutation probability P as a consequence of their quantum- 
mechanical nature. In superstring networks, intercommutation 
events are, therefore, less frequent, which demands longer 
numerical simulations (potentially up to a factor 10? longer). 
Such simulations have not yet been carried out, which makes 
the prediction of the GW signal from cosmic superstrings rather 
uncertain at present. In the absence of a more sophisticated 
treatment, we follow the standard practice in the literature and 
simply estimate the GW spectrum from cosmic superstrings by 
rescaling the spectrum in Equation (40), 

BD = 508 (52) 
where we allow for P values drawn from a log-uniform prior, 
logio Р Є [—3, 0], which covers the theoretically expected 
range for cosmic F- and D-superstrings (Jackson et al. 2005). 


5.4.2. Results and Discussion 


We now turn to the outcome of our Bayesian fit analysis, 
discussing first our results for stable strings. As evident from 
Figures 2 and 3, we find that stable strings are not able to 
provide a better fit of the PTA data than the benchmark SMBHB 
model. In fact, among all exotic GWB sources considered in 
this work, stable cosmic strings represent the only case that 
consistently yields a Bayes factor in favor of the SMBHB 
interpretation. Comparing the STABLE-C, STABLE-K, STABLE- 
M, and STABLE-N models to the SMBHB model, we obtain 
В = 0.277 = 0.006, В = 0.364 + 0.008, В = 0.379 + 
0.008, and В = 0.307 = 0.006, respectively; comparing the 
STABLE-C+SMBHB, STABLE-K+SMBHB, STABLE-M+SMBHB, 
and STABLE-N+SMBHB models to the SMBHB model, we 
obtain B = 0.76 + 0.01, В = 0.89 + 0.02, В = 0.84 + 0.02, 
and B = 0.83 + 0.01, respectively. These Bayes factors are 
close to unity, which means that adding a cosmic-string 
contribution to the GWB signal on top of the SMBHB signal 
does not improve the fit. Meanwhile, the larger prior volume 
pushes the Bayes factors to values slightly smaller than unity. 

The reason for the poor performance of the stable-string 
models is straightforward. In order to explain the relatively 
large amplitude of the signal, a comparatively large value of the 
cosmic-string tension Gj is necessary. Large Gy values, 
however, tend to result in a rather flat GW spectrum in the PTA 
band as seen in Figure 20 (see also, e.g., Figure 1 in Blanco- 
Pillado et al. 2011 or Figure 1 in Blasi et al. 2021), which 
clashes with the fact that the data seem to prefer a blue-tilted 
I? Ow spectrum. If we approximate the GW signal from stable 
strings by a simple power law, this observation can be 
reformulated in the language of the 7-А plot in Figure 1. In 
terms of y and A, we then conclude that stable strings allow one 
to obtain y < 4 only for log,, A < —15.0, while larger values in 
the range logo A ~ —(14.5 --- 14.0) can only be achieved for 
гу ^» 5. Тһе GW spectrum from stable strings is therefore always 
either too weak or too flat. 

In Figure 9, we present the reconstructed posterior distribu- 
tions for the cosmic-string tension in our four stable-string 
models. The left panel of Figure 9 shows our results in the 
absence of an additional SMBHB contribution to the signal, 
while the right panel displays the posterior distributions in the 
combined strings-plus-SMBHBs models. In the former case, 
we find peaked distributions centered around values of the 
order of logi) Gu ~ —(10.5 --: 10.0). Values of the cosmic- 
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Figure 9. Reconstructed 1D marginalized posterior distributions for the dimensionless cosmic-string tension parameter Си in the four stable-string models without 
SMBHBs (STABLE-C, STABLE-K, STABLE-M, and STABLE-N) in the left panel and including SMBHBs (STABLE-C+SMBHB, STABLE-K+SMBHB, STABLE-M+SMBHB, 
and STABLE-N--SMBHB) in the right panel. The dashed vertical lines show the 68% Bayesian credible intervals, which we construct as described in the text, and the 
solid vertical lines mark ће Си values where the К ratio in Equation (10) reaches a value of 1/10. Figure 28 in Appendix C.4 shows an extended version of this plot 


that includes the SMBHB parameters Авнв and Увнв: 


string tension in this range represent a compromise between the 
two extremes discussed above: at smaller log Gu the GW 
signal becomes too weak, while at larger log, Gu it becomes 
too flat. Moreover, we note that the order of the peaks in the 
posterior distributions—STABLE-M, STABLE-K, STABLE-N, 
STABLE-C from left to right—agrees with the ordering found 
in Blanco-Pillado et al. (2021) at y 4.6. The underlying 
reason is that, for fixed Си, STABLE-M predicts the strongest 
GW signal in the nanohertz frequency band, followed by 
STABLE-K and STABLE-N, while STABLE-C predicts the 
weakest signal. 

If we extend the stable-string models by including an 
SMBHB contribution to the GW signal, the posterior 
distributions feature not only a peak toward the upper end of 
the log; Gu range but also an extended plateau at small 
log; Си values. This plateau corresponds precisely to the 
plateau that we alluded to in the definition of the K ratio in 
Equation (10) and hence allows us to derive the following 
upper limits on the cosmic-string tension in the STABLE-C 
+SMBHB, STABLE-K+SMBHB, STABLE-M+SMBHB, and 
STABLE-N+SMBHB models, respectively: logo Gu < —9.67, 
—9.87, —10.10, and —9.71. We recall that, for cosmic-string 
tensions exceeding these values, the likelihood of the strings- 
plus-SMBHBs model drops below 1/10 of the likelihood of the 
SMBHB benchmark model. In this region of parameter space, 
adding GWs from stable strings to the signal does more harm 
than good and is strongly disfavored by the data. For 
comparison, we also quote the upper limits of the 95% 
Bayesian credible intervals for Gy, specifically of the 95% 
highest posterior density intervals (HPDIs). For each model, 
the HPDI is the narrowest possible range that includes 95% of 
the posterior probability. By construction, the boundaries of an 
HPDI are at points with the same posterior probability density, 
and each parameter value inside the HPDI has higher posterior 
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probability density than any parameter value outside the 
НРОТ.® We find log; Gu < —9.88, —10.04, —10.26, and 
—9.90 for STABLE-C+SMBHB, STABLE-K+SMBHB, STABLE-M 
+SMBHB, and STABLE-N--SMBHB, respectively. 

Next, we turn to metastable strings, which, as can be read off 
from Figure 2, provide a better fit to the data. In the absence of 
an SMBHB contribution to the GW signal, the Bayes factors 
for the META-L and META-LS models are B = 13.4 + 0.4 and 
B = 21.3 + 0.8, respectively. Adding an SMBHB contribution 
to the GW signal, the Bayes factors for the META-L and META- 
LS models are instead В = 11.1 + 0.3 and В = 18.9 + 0.7, 
respectively. Again, the fit does not become better, but the 
larger prior volume results in a decrease of the Bayes factors. 

The reconstructed 1D and 2D posterior distributions for the 
two metastable-string models are shown in the corner plots in 
Figure 10. From these plots, we read off that the NG15 
data prefer values of the decay parameter vV% of around 
«кос» 8 in combination with a large cosmic-string tension, 
log Gu ~ —(7 --- 4). Here Vk ~ 8 causes a suppression of 
the GW signal in the nanohertz band compared to stable strings 
because of the exponential factor in Equation (50). This 
suppression results in a decreased signal strength and a larger 
spectral tilt. The decrease in signal strength can, however, be 
compensated by a large cosmic-string tension, such that the 
final spectrum still has the right amplitude but is now more 
blue-tilted than in the stable-string case. This interplay of the 


36 Th the case of the four strings-plus-SMBHBs models, the probability density 
threshold for this procedure may coincide with the height of the plateau at low 
frequencies. If this happens, sampling fluctuations in the plateau region will 
lead to a set of disjoint intervals. To avoid these fictitious intervals, we set the 
upper limit where the posterior probability density falls to the level of the 
plateau and adjust the lower limit within the plateau so that the intervals contain 
95% of the posterior probability. 
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Figure 10. Same as in Figure 5 but for the metastable-string models META-L (GW emissions from loops only) in the left panel and META-LS (GW emission from loops 
and segments) in the right panel. The gray shaded regions are strongly disfavored by the NG15 data, as they result in a K ratio of less than 1/10 (see Equation (10); 
the teal shaded regions are ruled out by the CMB bound on cosmic strings (Ade et al. 2016; Charnock et al. 2016; Lizarraga et al. 2016), and the regions to the right of 
the yellow contour lines violate the LVK bound in Equation (24). Figure 29 in Appendix C.4 shows extended versions of the two plots that include the SMBHB 


parameters Авнв and Эвнв- 


different factors affecting the GW spectrum from metastable 
strings effectively leads to a better fit. 

Let us comment on a few characteristic features of the 
posterior distributions shown in Figure 10. First, we find that 
the log; Си posterior in the META-LS model extends to slightly 
larger values than in the META-L model, which demonstrates 
the effect of the additional contribution to the GW signal from 
string segments. In fact, in the region of highest posterior 
density, the segment contribution dominates over the loop 
contribution in the META-LS model. Second, we point out that 
the 1D marginalized posterior distributions for Іов Си exhibit 
small local maxima at log,, Gu ~ —(11 ес 10), which corre- 
spond to the stable-string limit within the metastable-string 
models. This limit is realized for large values of the decay 
parameter, VK > 9, which pushes the effect of the network 
decay to frequencies below the PTA band. Next, we observe 
that for META-L the logio Си posterior experiences a sharp 
drop-off at log, Си ~ —5, whereas for META-LS there is а 
small dip in the log, Gu posterior at logi Gu ~ —5. Both 
features can be traced back to the Heaviside theta function in 
Equation (46), which ensures that no more new loops are 
formed during the decay regime of the network. Because of this 
Heaviside theta function, the loop contribution to the GW 
spectrum moves to frequencies above the PTA band if we raise 
log; Gu above log, Gu ~ —5. In this sense, the drop-off and 
the dip in the log,, Си posteriors should be regarded as being 
due to our simplified theoretical modeling of the GW spectrum. 
More work is necessary to improve on the description in terms 
of a simple Heaviside theta function and obtain a better 
understanding of the evolution of the decaying network. We 
expect that a more accurate description of the transition from 
the scaling to the decay regime would result in smoother 
log; Си posteriors. Finally, we note that some of the 
fluctuations in the posteriors in Figure 10 can be attributed to 
the fact that, in the case of the metastable-string models, we 
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have to work with tabulated data for the GW spectrum, based 
on the numerical evaluation of Equation (40). 

We also use the corner plots in Figure 10 to highlight the 
relevant bounds on the parameter space spanned by logi) Си 
and vK. The gray shaded areas notably indicate the K-ratio 
bound, which marks the parameter regions that are ruled out by 
the NANOGrav data. In these regions, the GW signal from 
metastable strings exceeds the observed signal and hence is 
unacceptably large. We find that the NANOGrav bound is 
stronger than the well-known CMB bound, which demands that 
a cosmic-string network that has not yet decayed by the time of 
recombination must not have a cosmic-string tension larger 
than log; Gu = —7 (Ade et al. 2016; Charnock et al. 2016; 
Lizarraga et al. 2016). In order to derive the CMB bound, we 
estimate that a decaying network completely disappears 
because of GW emission at X times around 
t, c (2/(ГСиу)! t, which is the time when the second term 
in the exponent in Equation (50) becomes large. The teal 
shaded areas in Figure 10 indicate where the conditions 
log Gu > —7 and te > tres = 370,000 yr are satisfied simulta- 
neously. Last but not least, the yellow solid lines represent the 
LVK bound in Equation (24), which translates to the upper 
limit Си <2 x 10 7.77 

The LVK bound on the amplitude of the stochastic GWB 
appears to be in tension with most of the parameter space 
preferred by the NANOGrav data. In particular, the 6896 
credible regions in the log) Gu—/K parameter plane are mostly 


37 In our analysis, we work with the upper limit on the amplitude of a generic 
isotropic and stochastic GWB that was reported by the LVK Collaboration in 
Abbott et al. (20214). Combining this limit with our own GW spectra, we find 
Си < 2 x 10 7. This bound needs to be compared to the results of Abbott 
et al. (2021b), a dedicated search for GWs from stable cosmic strings by the 
LVK Collaboration. The analysis in Abbott et al. (2021b) models the GW 
spectrum from cosmic strings based on slightly different assumptions. Model A 
in Abbott et al. (2021b) is, however, similar to our approach and leads the same 
bound on Си in the limit of a large number of kinks, № > 1 (see the panel for 
model A in Figure 3 of Abbott et al. 2021b at №, ~ 100 --- 200). 
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ruled out by this bound. However, we stress that the LVK 
bound in Figure 10 relies on an extrapolation of the GW signal 
across 10 orders of magnitude in frequency, from РТА 
frequencies, f~ few х 10 ?nHz, to LVK frequencies, 
f^ few x 10 Hz. In order to perform this extrapolation, we 
have to assume a cosmological expansion history across 10 
orders of magnitude in temperature, even though not much is 
known about the equation of state of the universe prior to BBN. 
Any deviation from the expansion history of standard big bang 
cosmology can therefore affect the extrapolation of the GW 
spectrum to higher frequencies and potentially render the LVK 
bound harmless. A simple example is an early stage of matter 
domination (Allahverdi et al. 2021), which would suppress the 
GW signal from metastable strings at high frequencies and thus 
allow for the possibility of cosmic-string tensions as large as 
those in Figure 10. On the other hand, if we trust the 
extrapolation to higher frequencies, we conclude that some 
parts of the 95% credible regions in Figure 10 are in fact not 
ruled out. In this case, metastable strings with «/к ~ 8 and 
log; Gu ~ —7 are still able to provide a good fit of the data, 
which represents a particularly interesting scenario for two 
reasons: a value of logo Gu ~ —7 would point to a GUT 
origin of the string network, and ground-based interferometers 
would be poised to detect the stochastic GWB signal from such 
a network in the near future. 

Finally, we discuss our results for cosmic superstrings, for 
which we obtain the largest Bayes factors among all cosmic- 
string models considered in this work: В = 46 = 2 for GWs 
from superstrings alone and В = 37 + 2 for an SMBHB 
contribution added to the superstring GW signal, both 
compared to the SMBHB model. 

Superstrings result in a good fit of the data because, unlike 
ordinary stable strings, they permit a GW signal that is neither 
too weak nor too flat. To see this, note that small cosmic-string 
tensions, log,, Gu < —11, yield a blue-tilted №Осу spectrum 
at nanohertz frequencies (see our discussion above). In the case 
of ordinary strings, the amplitude of this spectrum would be too 
weak; however, for superstrings the suppression of the GW 
signal at log; Gu 5 —11 can be compensated by the 1/P 
enhancement factor in Equation (52). By choosing P appro- 
priately, the amplitude of the GW signal can then be raised 
until a good fit of the data is reached. This interplay between 
the two parameters of the SUPER model can also be observed in 
Figure 11. The 2D posterior distribution for log,, Gu and P in 
this corner plot displays a strong covariance in line with our 
heuristic understanding. 

The highest posterior density is achieved at small intercommu- 
tation probabilities and cosmic-string tensions, log, P ~ —3 and 
logo Си ~ —12, where log, Р = —3 corresponds in fact to the 
edge of the prior range for P. This parameter region is in tension 
with the LVK bound in Equation (24) (see the yellow solid line in 
Figure 11) but may be viable assuming a nonstandard expansion 
history at high temperatures. Similarly to what is found in 
Figure 10, we highlight the region ruled out by the NGI15 data 
because it leads to a K ratio of less than 1/10. As expected, this 
region is located to the lower right of the 68% and 95% credible 
regions, where small log, Р values cause the amplitude of the 
GW spectrum to exceed the measured strength of the signal. The 
tension of cosmic superstrings can also be constrained by CMB 
observations (Charnock et al. 2016). Existing upper limits are, 
however, located at larger Си values and thus not relevant for our 
results in Figure 11. 
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Figure 11. Same as in Figure 10, but for the SUPER model. Figure 30 in 


Appendix C.4 shows an extended version of this plot that includes the SMBHB 
parameters Авнв and увнв. 


In conclusion, we stress again that the GW spectrum from 
cosmic superstrings is not well understood at present, and more 
work is needed to arrive at a reliable prediction. Such an effort 
will be important to confirm that GWs from cosmic super- 
strings do indeed represent a realistic and viable interpretation 
of the РТА signal. 


5.5. Domain Walls 
5.5.1. Model Description 


Domain walls are 2D topological defects that form when a 
cosmological phase transition results in the spontaneous 
breaking of a discrete symmetry, filling the universe with 
patches in different degenerate vacua(Kibble 1976). In the 
absence of significant interactions with relativistic particles, 
domain wall networks are expected to reach a scaling regime in 
which each Hubble volume, H^, contains A ~ 0(1) domain 
walls (see, e.g., Press et al. 1989; Hindmarsh 1996, 2003; 
Garagounis & Hindmarsh 2003). In this regime, the energy 
density stored in domain walls is given by 


Оруу ~ A cH, (53) 


where A is nearly constant and с is the domain wall tension, 
which gives the domain walls their surface energy density. 
During the scaling regime, the energy density stored in domain 
walls dilutes more slowly than that of relativistic radiation and 
nonrelativistic matter. Indeed, the total energy density of the 
background always scales like р. ос H°, which means that 
Оруу = ppw/ pe grows like the Hubble radius when the domain 
wall network is in the scaling regime, Qpw с H '. Therefore, 
domain walls eventually overclose the universe and alter the 
cosmological evolution in a way that is incompatible with 
CMB observations (Zeldovich et al. 1974). A possible solution 
to this problem is to have domain walls decay at some 
temperature T by assuming, e.g., that the global symmetry 
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responsible for domain wall formation is explicitly broken. 
Indeed, an explicit breaking of the symmetry introduces a bias 
among the possible low-energy vacua, thus lifting their 
degeneracy, which eventually leads to the collapse of the 
domain wall network (Vilenkin 1981b; Gelmini et al. 1989; 
Larsson et al. 1997). 

During the scaling regime, domain walls continuously 
change their configuration and shrink in order to maintain the 
scaling relation in Equation (53). While these processes take 
place, a fraction of the domain wall energy is released in the 
form of GWs, which produce a GWB with a present-day relic 
abundance given by(Vilenkin 1981b; Preskill et al. 1991; 
Gleiser & Roberts 1998; Chang et al. 1999; Hiramatsu et al. 
2010; Kawasaki & Saikawa 2011) 


WOGw(f) = ==? гаї S/P), 


T 


(54) 


where D is the dilution factor defined in Equation (35), € — 0.7 
is an efficiency coefficient derived from numerical 
simulations (Hiramatsu et al. 2014), and a,, is the fraction of 
the total energy density stored in domain walls at Ту, 


Ppw | 
зн?М |. 


ax = Ору (Тю) = (55) 


The GWB is dominated by the emission taking place just 
before the decay of the domain wall network. As the GWs are 
predominantly sourced by horizon-scale structures in the 
network, the typical frequency of the GWs emitted around 
this time is Ну, i.e., the Hubble scale at Tą. Therefore, after 
redshifting until today, the GWB is expected to exhibit a peak 
frequency 


10.75 d 
f, = 1.14 nHz| — | 
Ex,s 


8% )( Tx 
10.75 10 MeV 


| (56) 


with both gą and g.., evaluated at Tą. For the spectral shape, 
S(x), we use a parameterization similar to the one used for the 
phase transition spectrum in Section 5.3:5* 
(a + Б)" 
(bx-4/¢ ES ах?/‹)с 


S(x) = (57) 
Causality fixes а = 3, while numerical simulations (Hiramatsu 
et al. 2014) suggest b — c = 1. In our analysis, we therefore fix 
the low-frequency slope of the domain wall signal to the value 
predicted by causality, and, following Ferreira et al. (2022), 
we allow b and c to float within the uniform prior ranges 
b € [0.5, 1] and c € [0.3, 3]. 

In our analysis we consider two possible decay channels for 
the domain wall network: dark radiation (DW-DR) and SM 
particles (DW-SM). In the DW-DR model, instead of using o, 
we parameterize the strength of the domain wall signal in terms 
of the dark radiation produced in the decay, ров. This quantity 
is usually expressed in terms of the effective number of extra 
neutrino species, ANG, which we defined in Section 5.1 and 
which is bounded from above by BBN and CMB observations, 
ANG ~ few x 0.1 (Pisanti et al. 2021; Yeh et al. 2021). In 


88 Since the efficiency coefficient ё is derived by matching the value of the 
domain wall spectrum at the peak frequency, Equation (57) does not contain 
the normalization coefficient M appearing in Equation (38). 
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our MCMC analysis of the DW-DR model, we follow Ferreira 
et al. (2022) and use АМ.” = 0.39 as the upper bound on our 
prior for the parameter ДМ, ғ. Assuming that all the domain 
wall energy is converted to dark radiation at Tą, AN. is 


related to a, Бу 


1/3 
А№ = | Em | Gal Е ) 
Sx 8% 0.1 


where, as before, both g, and 8 ; are evaluated at Ту. 

In the DW-SM model, BBN restricts the possible values of the 
decay temperature to T, 2 2.7 MeV (Jedamzik 2006; Bai & 
Korwar 2022) for any detectable value of o. Following 
Ferreira et al. (2022), we also impose a, < 0.3 to avoid any 
possible deviation from radiation domination and to evade 
bounds from ANGE. 


(58) 


5.5.2. Results and Discussion 


The reconstructed posterior distributions for the parameters 
T, and o, (Ty, and ДМ) of the DW-SM (DW-DR) model are 
reported in Figure 12, for both the case where the domain walls 
are assumed to be the only source of GWs (blue contours) and 
the scenario where we consider the superposition of the domain 
wall and SMBHB signals (red contours). Full corner plots 
including the posterior distributions of the spectral shape 
parameters b and c are reported in Figure 31 in Appendix C.5. 

For both the DW-SM and DW-DR models, with and without 
the inclusion of the SMBHB signal, we find that the GWB 
produced by domain walls peaks around 10-8 Hz such that 
most of the low frequency bins are fit by the low-frequency tail 
of the spectrum (see Figures 3 and 19). Specifically, for the 
DW-SM (DW-DR) model, we find that 7, € [110, 275] ([79, 
153]) MeV at the 68% credible level and Ту Є [76, 505] ([54, 
198] MeV at the 9596 credible level. When including the 
SMBHB contribution to the GWB, we find T, € [108, 309] 
([54, 216]) MeV at the 68% credible level and Ту Є [67, 843] 
(no bound on T4) MeV at the 95% credible level for the DW-SM 
(DW-DR) model. We notice that, with and without the inclusion 
of SMBHBs, the recovered transition temperature for the DW- 
SM model is high enough to evade BBN constraints. 

For the DW-DR model, the posterior distribution for AN, is 
peaked near the upper prior boundary, signaling that larger 
values fit the observed signal better. Specifically, we find 
AN. > 0.32 at the 68% credible level and AN, Z 0.25 at the 
95% credible level. Including the contribution from SMBHBs 
allows the distribution for AN, to extend to lower values, and 
we find AN, > 0.23 at the 68% credible level and no bound at 
the 95% credible level. We thus conclude that the DW-DR 
model prefers large ANG values in the vicinity of existing 
bounds. This means that the most promising parameter regions, 
i.e., regions that аге not yet ruled by AN. but still manage to 
fit the NANOGrav signal, point to ANG values within the 
reach of upcoming experiments, including CMB-S4 (Abazajian 
et al. 2022), which promises to be sensitive to AN. values as 
small as AN, 0.06 at the 95% confidence level. 

For the DW-SM model, we find а, Є [0.080, 0.19] at the 68% 
credible level and сих € [0.053, 0.35] at the 95% credible level. 
With the inclusion of the SMBHB contribution, smaller values 
of аж become allowed, and we find o Є [0.079, 0.22] at the 
68% credible level and œx € [0.047, 0.61] at the 95% credible 
level. We also notice a partial degeneracy between Tą and Qx, 
due to low frequency bins (which are the ones contributing the 
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Figure 12. Same as in Figure 5, but for domain walls decaying to SM particles (left panel) and a dark sector (right panel). Figure 31 in Appendix C.5 shows extended 
versions of the two plots that include the spectral shape parameters b, c and the SMBHB parameters Авнв and увнв. 


most to the evidence for a GWB) being fit by the low-frequency 
tail of the domain wall spectrum. Therefore, a shift to higher 
transition temperatures (and therefore to higher peak frequen- 
cies) can be partially reabsorbed with a shift to higher o. values. 

For both DW models, we identify regions of the parameter 
space where the GW signal from domain walls is too strong to 
be compatible with NG15 data. These regions, shaded in gray 
in Figure 12, illustrate for the first time how PTA data can be 
used to constrain models of domain walls. 

Finally, we state the Bayes factors for the model comparison 
between the domain wall models and the SMBHB reference 
model. The Bayes factors for DW-SM and DW-DR versus 
SMBHB are В = 14.8 + 0.5 and В = 1.62 + 0.05, respec- 
tively, while the Bayes factors for DW-SM--SMBHB and DW-DR 
-SMBHB versus SMBHB are 6B=21.1+0.9 апа 
В = 2.53 + 0.10, respectively. In both cases, the extra 
SMBHB contribution helps to improve the fit of the NGI5 
data. In the case of DW-DR, we moreover observe the same 
effect as for PT-SOUND in Section 5.3: adding the SMBHB 
contribution to the GWB signal results in a plateau region in 
the posterior distribution of the DW-DR model parameters that 
manifests itself as part of the 95% credible region in Figure 12. 


6. Deterministic Signals from New Physics 


In addition to the GWB signals discussed previously, there 
are several new-physics theories that can imprint a determinis- 
tic signal, described by a time series h, in pulsar timing data. In 
this section, we consider the deterministic signals induced by 
ULDM and DM substructures. After finding no statistically 
significant evidence for such signals in our data, we report 
upper limits on the allowed strength of these signals. 


6.1. Ultralight Dark Matter 


6.1.1. Model Description 


While DM constitutes roughly 27% (Zyla et al. 2020) of the 
energy density of the universe, very little is known about its 
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fundamental properties. Consequently, a wide range of DM 
models remain consistent with cosmological and astrophysical 
observations. The lightest possible DM particles are classified 
as ultralight, or fuzzy, DM. These particles must be bosonic; 
otherwise, they could not be packed into galaxies owing to 
Fermi degeneracy pressure (Tremaine & Gunn 1979; Di Paolo 
et al. 2018; Savchenko & Rudakovskyi 2019; Alvey et al. 
2020; Davoudiasl et al. 2021). These ULDM models also 
generically suppress structure on small scales, allowing them to 
be potential solutions to the small-scale structure problems of 
the standard ACDM paradigm (Bullock & Boylan-Kol- 
chin 2017). However, too much suppression on large scales 
would be in conflict with CMB measurements (Hlozek et al. 
2015; HloZek et al. 2017, 2018), which sets a lower bound on 
the DM mass, 107^ eV < mo. PTAs can probe these miniscule 
masses, since, as we describe below, the frequency of the 
ULDM signal, f, is generally proportional to the ULDM mass, 
2mf ~ mg. Therefore, the sensitivity window of NANOGrav is 
expected to be 1072? еу < ть 5 10720 eV. 

Other astrophysical constraints, such as measurements of ће 
Гуа forest (Armengaud et al. 2017; Iršič et al. 2017; Kobayashi 
et al. 2017; Rogers & Peiris 2021), galactic subhalo mass 
functions (Schutz 2020; Banik et al. 2021; Nadler et al. 2021), 
and stellar kinematics (Dalal & Kravtsov 2022), can push this 
bound further up, ranging from 1072! to 10 P eV, and the 
nonobservation of superradiance at SMBHs (Arvanitaki et al. 
2015; Stott 2020; Ünal et al. 2021) pushes the bound up to 
1077! ...107" eV. РТА searches for ULDM іп the 
10723 еу < т, X10 ? eV window can then be viewed as 
complementary to these astrophysical searches. Signals in 
PTAs do not depend on the same astrophysical uncertainties, 
e.g., modeling the nonlinear small-scale matter power spectrum 
with analytic or numeric methods (Zhang et al. 2018, 2019) or 
the evolution of specific density profiles. Moreover, the power 
of these methods quickly degrades when considering sub- 
components of DM, or populations of DM that do not compose 
the whole DM density. The searches discussed here have a 
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Table 1 
Summary of the Different Effects Induced by ULDM That Generate a Deterministic Signal of the Form Given in Equation (59) 
Effect ALG) Ab у(х) w = 2nf Spin (55) 
A A тарь n "Gp, na 

Metric fluctuations (28, +1) 2.3 Pe) QS, + 1) x Q^) 2m, 0,1 

— 205 a A [Po ^ ^ А 
Doppler UC) forces 10 ` e) Q(x) P Me nd (пу - є) ф(х) то l 
Pulsar spin fluctuations 0 ypdp ү = (x) mo 0 
Reference clock shift A 0 m 0 

Ypde m2 ф(х) Ф 


Note. “Е” and “Р” subscripts denote Earth and pulsar term contributions, respectively. A‘ is the signal amplitude for ће ith DM polarization, wis ће signal frequency, 
and the "Spin" column refers to whether the effect occurs for scalar (spin-0) or vector (spin-1) ULDM candidates. p, is the local DM density, taken to be 
0.4 GeV cm ^, ть is the ULDM mass, g are gauge couplings, О is the charge under the corresponding gauge group, y parameterizes the sensitivity to different 


fundamental constant fluctuations, d are defined in Equation (64), and A = Mpi/ Ап. ф(х) is a random variable representing the fluctuations in the ULDM field in 
different correlation patches; the correlated limit corresponds to ф(х) — Ф. n, is a vector pointing from Earth to the Аһ pulsar, and є; аге the ULDM polarization 


vectors. 


weaker dependence on the subcomponent fraction and are 
therefore still useful in the hunt for DM. 

While there is a large phenomenology of ULDM signals that 
can be produced in PTAs, the timing residuals for the /th 
pulsar, Л, are deterministic and can be written in the form 


тб) = 5 7 Ag Gn) sinat + УЕ) 


+ Ab (xp,)sin(wt + уБ), (59) 


where i indexes the DM field polarization, "? w is the signal 
frequency, and we have split the signal into two terms: an 
"Earth term" with amplitude Ag (xg) and time-independent 


phase Y, and a “pulsar term” with amplitude Abi (хр) and 


time-independent phase т г The phases can be written as 
(60) 
(61) 


Ve = arg), 
Урі = «хе — xp | + o (Хр), 
where о, '(х) is dependent on the underlying ULDM field phase 


at point x, and the additional term in Ve 1 15 due to the light- 


travel time between Earth and the pulsar. 
An important scale in understanding these signals is given by 
the correlation length of the ULDM field, £e, 


10-22 eV | 


то 


2m 


ov» 


ES ~ 0.4 d (62) 


where v~ 10 ? is the ULDM velocity. If the correlation 
length is much larger than the distance between Earth and the 
pulsar, £. >> [хе — хь, both experience the same DM field. In 
this "correlated" limit of the signals, the amplitudes and phases 
can be taken to be position independent: 


Ар (хв) — Ар, 


Ap(xp;) > Арі, а) > а. (63) 


Generally, a correlated analysis can drastically reduce the 
number of free amplitude parameters in the MCMC. For 
example, if ULDM couples to all pulsars identically, then there 


89 Here we focus on signals from scalar (зрш-0) and vector (spin-1) ULDM. 
We leave the search for the effects of spin-2 ULDM (Marzola et al. 2018; 
Armaleo et al. 2020a, 2020b; Xia et al. 2023) to future work. For scalar DM, i 
is absent, since the field has only one component. For vector DM, i € (1, 2, 3} 
for each massive vector polarization. 
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is only one amplitude, Ap, instead of one for each pulsar, to 
account for fluctuations in the ULDM field. However, since 
тхе — Xp; 2 1, ye, and ^g must still be taken to be 
independent. In the uncorrelated limit, 4. «Є [хь — xp,|, the DM 
field is no longer correlated, and its amplitude fluctuations in 
different patches need to be accounted for. 

In Table 1, we summarize the ULDM signals searched for in 
our analysis by providing the corresponding expression for the 
amplitudes and signal frequency appearing in Equation (59). In 
the following, we give a brief review of each of these signals 
and refer the reader to the original references for more details: 


1. Metric fluctuations (Khmelnitsky & Rubakov 2014; Por- 
ayko & Postnov 2014; Porayko et al. 2018; Kato & 
Soda 2020; Nomura et al. 2020; Unal et al. 2022; Wu 
et al. 2022): Oscillations in the ULDM field generate 
fluctuations in the local stress-energy tensor, 7,,,, which 
are independent of any direct couplings the ULDM may 
have with SM fields. These fluctuations in 7,,, generate 
fluctuations in metric perturbations Бу  Einstein's 
equations. These metric perturbations then affect the 
photon geodesic on its path from the pulsar and generate 
a timing residual with pulsar and Earth amplitudes given 
in Table 1. While the amplitudes shown in Table 1 are 
specific to scalar (spin-0) ULDM, vector (spin-1) ULDM 
also generates this purely gravitational signal. The vector 
DM scenario is more complicated (Nomura et al. 2020), 
as off-diagonal components of the metric fluctuations can 
generate nontrivial angular correlations between signals 
seen in different pulsars. In our analysis, we ignore these 
spatial correlations and set approximate bounds by 
treating the vector as three scalar components, such that 
the vector signal is three times larger than the scalar 
one (Nomura et al. 2020). 

2. Doppler-U(1) forces (Graham et al. 2016; Xue et al. 
2022): A background of ultralight vector DM can create 
dark "electric" fields, which can accelerate pulsars, or 
Earth, and Doppler-shift light arrival times. For example, 
if the DM is the gauge field of a local baryon symmetry, 
U(1)s, then it couples to the Earth and pulsar baryon 
current density. This coupling generates a force on both 
Earth and pulsars proportional to their baryon number, 
Ов = № х qp, where Ng is the number of baryons and дв 
is the gauge charge. This scenario is straightforwardly 
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generalized to other U(1) models, e.g., models where the 
difference between baryon and lepton numbers, B — L, is 
a local symmetry. These forces cause periodic displace- 
ments between Earth and pulsars and generate timing 
residuals with amplitudes given in Table 1. These forces 
also exist for scalar DM coupling to SM operators (e.g., 
р fin, where ф is the DM and п is the neutron field). 
However, in the scalar case, these forces originate from 
the field gradient Уф (Graham et al. 2016), which causes 
them to be velocity suppressed compared to the vector 
DM scenario. 


The next two signals we discuss are specific to scalar 
ULDM. All the linear couplings of a scalar ULDM field to the 
SM can be summarized in a single Lagrangian, 


d d 
Fy FY + «03 
4е 283 


£3 


x Gi, Gi" 


- У) ат; У) (а + чадтаа|, | (64) 


1=е, и 4=и,ӣ 


where Е,,, and G,,, are the photon and gluon field strengths, 
respectively, Л = Мы/ Ап, 33 is the QCD beta function, ^fa 
are the light quark anomalous dimensions, ту are the fermion 
masses, and the d values are dimensionless couplings. These 
couplings induce periodic oscillations in the values of 
fundamental constants, i.e., particle masses and couplings, 
which can affect timing residuals in two ways: 


1. Pulsar spin fluctuations (Kaplan et al. 2022): Particle 
mass fluctuations change the moment of inertia of the 
pulsar, which, by conservation of angular momentum, 
leads to pulsar spin fluctuations. The amplitude of the 
signal is given in Table 1. The y parameters are pulsar 
specific and denote the pulsar sensitivity to a specific 
coupling, e.g., y, parameterizes the pulsar sensitivity to 
changes in the electron mass. у„ parameterizes the pulsar 
sensitivity to the mass-weighted combination of quark 
mass couplings, ад = (т.а, + mgdg)/(m, + та). To 
be explicit, following Kaplan et al. (2022), we employ in 
our analysis 


1071, 
1075, 


#=-5, 


у. = 2 x 103, 


ув = —24 x 


y217 x (65) 


for each pulsar, which assumes the simplest possible 
model for the pulsars' moment of inertia. 

2. Reference clock shifts (Graham et al. 2016; Kaplan et al. 
2022): PTA timing residuals are measured with respect to 
a collection of mostly cesium (McCarthy 2009) atomic 
clocks. Therefore, changes to these atomic clock 
frequencies can appear in PTA data as apparent shifts 
in timing residuals. These shifts are perfectly correlated 
among pulsars and have the amplitude reported in 
Table 1. Similar to the pulsar spin fluctuation signal, 
the y parameters denote the atomic clock sensitivity to 
specific couplings. To be explicit, following Kaplan et al. 
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(2022), we take these parameters to be 
=1 ya = 0.158 
y,=2 y, = 4.83, (66) 


assuming cesium clocks to be used as a reference. 


6.1.2. Results and Discussion 


We search for ULDM signals by performing a variety of 
Bayesian fit analyses on the NGI5 data set. The priors for the 
ULDM model parameters are summarized in Table 3, while the 
priors for the intrinsic red-noise parameters are reported in 
Table 2. Since in this analysis we want to remain agnostic 
about the origin of the GWB, we model the GWB as a power 
law and allow the values of the amplitude and spectral index to 
float within the following prior ranges: logjjAGws Є 
[—18, —11] and *ows Є [0, 7]. 

Over most of the ULDM mass range, 10 "eV < 
mg 5 10 7? eV, we find no significant evidence for any of 
the ULDM signals described in the previous section (with 
B ~ 1 іп favor of models including a ULDM signal on top of 
the GWB). However, for ULDM models that give rise to an 
"Barth term" signal, we find mild evidence for an ULDM signal 
with frequency f-4nHz. Specifically, restricting the prior 
range to f € [3.05, 6.09] nHz, we find a Bayes factor of B ~ 2 
(B ~ 1.5) in favor of the model including the ULDM signal in 
the correlated (uncorrelated) regime. This result is expected 
given the excess of monopolar-correlated power observed in 
the second frequency bin of the common red-noise process (see 
the discussion in NG15gwb for more details). Unfortunately, an 
ULDM interpretation of this monopolar signal is difficult, since 
the corresponding ULDM masses, m, ~ 2 х 10723 eV, needed 
to explain this excess are in tension with other astrophysical 
bounds: Гуа forest (Armengaud et al. 2017; Iršič et al. 2017; 
Kobayashi et al. 2017; Rogers & Peiris 2021), galactic subhalo 
mass functions (Schutz 2020; Banik et al. 2021; Nadler et al. 
2021), and stellar kinematics (Dalal & Kravtsov 2022). 

Without convincing evidence for a signal, we compute 
constraints оп the ULDM model parameters, shown in 
Figures 13-15. All constraints are the 95th percentile of the 
marginalized posterior distribution for the parameter on the 
vertical axis. The curves labeled “(un)correlated” correspond to 
the analysis done in the (un)correlated limit, discussed in the 
previous section. 

In Figure 13, we show the constraints on the local ULDM 
energy density that can be derived assuming only gravitational 
coupling between the ULDM and SM fields. In Figure 13, we 
show the constraints on the local ULDM energy density that 
can be derived assuming only gravitational coupling between 
the ULDM and SM fields. The strongest bounds are obtained in 
the mass range mg < 1072? eV, where we nearly constrain 
ULDM to be a subcomponent of the total DM abundance. 
While we show constraints down to m, = 107^ eV, it is easy 
to extrapolate them to lower masses, where we expect them to 
remain flat down to т ~ 10 7 eV, where 1 /m, becomes of the 
same order of the interpulsar separation and the ULDM signal 
is additionally suppressed (Khmelnitsky & Rubakov 2014; 
Unal et al. 2022). While future PTA analyses will be able to 
improve on these constraints, constraining ULDM with 
тд 2, 1072? еу to have a local abundance smaller than 
рь = 0.4 GeV cm ? will be challenging, given that the 


constraints scale as pin x m3 for ть Z 1/Tobs. 
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Figure 13. Constraints on the local ULDM density for the correlated (solid 
line) and uncorrelated (dashed line) signals at the 9596 credible level (see the 
discussion in the text) The gray dashed line indicates the predicted DM 
abundance. 
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In Figure 14, we show the constraints for all the а parameters 
describing the scalar ULDM Lagrangian in Equation (64). Note 
that, for each panel, we assume that the parameter constrained 
is the only nonvanishing one. Overall, the limits are in rough 
agreement with the projections from Kaplan et al. (2022) and 
competitive with laboratory constraints (Boulder Atomic Clock 
Optical Network et al. 2005; Hees et al. 2016; Bergé et al. 
2018; Kumar Poddar et al. 2019; Dror et al. 2020). The 
strongest constraints, relative to laboratory bounds, are for 
ULDM models coupled to the electron, de, and muon, d,,, mass 
terms. Indeed, relative shifts of the energy levels utilized in 
atomic clock experiments are insensitive to the d, coupling, 
since atomic energy levels in different atoms scale identically 
with electron mass, leading to no relative energy level 
shifts (Van Tilburg et al. 2015; Kaplan et al. 2022). Such an 
insensitivity is not a problem for the PTA observable, though, 
since the pulsar phase evolution is not affected in the same way 
as the atomic energy levels. The lack of laboratory constraints 
оп d, is simply because there are not a large number of muons 
to study on Earth, whereas pulsars host a large number of them. 
As for the gravitational signal, we can extrapolate the 
constraints to lower masses, where we expect them to scale as 
ах 1/т down to m;-10 ?6eV, where 1/m, becomes 
comparable to the interpulsar separation. 

In Figure 15, we show constraints on the gauge coupling of 
models where Ше ULDM is the gauge boson of either U(1)g or 
СО )вог. Our constraints are roughly consistent with those 
published by the PPTA Collaboration (Porayko et al. 2018). 
This result is somewhat expected: while NANOGrav observes 
more pulsars, the average observation time is longer in PPTA, 
so roughly similar bounds are expected. 

The constraints presented in Figures 14 and 15 assume that 
ULDM makes up the entire DM content of the universe. 
However, if ULDM is only a subcomponent of the total DM 
abundance, these constraints can be easily rescaled. Indeed, 
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from Table 1, we see that the amplitudes for the direct coupling 
signals scale as d /р, for the scalar case апа О /р, for the 
vector case. Therefore, if ULDM is only a faction f of the total 
DM abundance, the constraints are weakened by a factor > 7 


Lastly, we comment оп the prior choice for Ф. Previous 
studies (Porayko & Postnov 2014; Kato & Soda 2020) assumed 
that ф is not a random variable in the problem and placed 
constraints assuming that this parameter is simply 1. However, 
as pointed out by Centers et al. (2021), this assumption is not 
appropriate. Since the observation time of PTAs is much 
smaller than the coherence time of the DM, т ~ (mv?) 1, only 
one instance of the field is being sampled within each 
correlation patch. The DM density in any correlation patch is 
then a random number, which follows a Rayleigh distribution 
with mean pọ, and the priors should be chosen to reflect this. 
We also note that while Porayko et al. (2018) find the 
distribution of through numerical simulation, these results are 
consistent with the analytic predictions by Foster et al. (2018). 


6.2. Dark Matter Substructures 
6.2.1. Model Description 


In the ACDM model, the structures we observe in the 
universe are seeded by primordial curvature fluctuations 
generated during inflation and then imprinted onto the DM 
density field. CMB observations indicate that these fluctuations 
have a nearly scale-invariant power spectrum on large scales ( 
i.e., for comoving wavenumbers kc Mpc '). However, on 
smaller scales, various theories of DM leave unique fingerprints 
on primordial perturbations or their evolution, resulting in 
different predictions for the amount of subgalactic DM 
substructures. Consequently, measuring local DM substruc- 
tures could be crucial in determining the correct model of DM. 

PBHs are perhaps the simplest example of such small-scale 
DM substructures. They can be formed in inflationary theories 
that create large density fluctuations on small scales (like the 
ones described in Section 5.2). Several studies investigated the 
possibility of identifying a galactic PBH population by 
analyzing the Doppler and Shapiro signals they can leave in 
PTA data (Seto & Cooray 2007; Siegel et al. 2007; Dror et al. 
2019; Ramani et al. 2020; Lee et al. 2021a, 2021b). In this 
analysis, we will closely follow the method outlined by Lee 
et al. (2021b) to constrain the local PBH abundance.” 

The Doppler signal results from the apparent shift in the 
pulsar spin frequency, generated by the acceleration induced by 
the gravitational pull of a passing PBH. According to the 
timescale of the transit event, 7, the signal can be further 
classified as dynamic (static) if T is much smaller (larger) than 
the observation time, Tops. In the static regime, the leading- 
order term of the Doppler signal that is not degenerate with the 
timing model is given by (Ramani et al. 2020; Lee et al. 
2021a, 2021b) 


Abs 
„sta 13, 


hp, sta (t) = 2 
yr 


(67) 


where Ap sta is a dimensionless parameter that can be related to 
the kinematic parameters of the transiting event (see 


90 А similar approach could be applied to set constraints on the local 
abundance of DM subhalos. However, we do not consider this case, since our 
constraints for PBHs are already quite weak. Constraints on DM subhalos 
would likely be even weaker, making it a less promising target for future PTAs. 
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Figure 14. The red solid (dashed) lines show our constraints at the 9596 credible level on the d parameters of the ULDM model, Equation (64), derived by searching 
for an (un)correlated signal. We assume that all reference clocks use the transition between the two hyperfine levels of the Св ground state. Additionally, for each 
panel we assume that the d parameters not shown are set to zero. The lower gray shaded regions correspond to regions of parameter space where the signal amplitude 
is less than the purely gravitational signal. Current constraints *Rb/Cs atomic clocks" (purple) are from Hees et al. (2016), “Al/Hg atomic clocks" (turquoise) are from 
Boulder Atomic Clock Optical Network et al. (2005), “MICROSCOPE” (teal) are from Bergé et al. (2018), “H/Si clock shift" (orange) are from Boulder Atomic 
Clock Optical Network et al. (2005), and “NS binary system" are from Kumar Poddar et al. (2019) and Dror et al. (2020). 


Appendix C.6 for more details). In the dynamic limit, and 
assuming that the signal is dominated by the closest transiting 
PBH, we get 


hp,ayn(t) = Ap,dyn (t = to) Ot Б 10), (68) 


where © is the Heaviside step function, Ара, is а 
dimensionless amplitude that can be related to kinematic 
parameters of the transiting event, and ig is the time of closest 
approach (see Appendix C.6 for more details). 

The Shapiro signal refers to shifts in the TOAs caused by 
metric perturbations along the photon geodesic produced by 
РВН$ transiting along the observer's line of sight. In the static 
limit, and after subtracting away degeneracies with timing 
model parameters, the leading-order term of this signal can be 
parameterized as 


As, t 
hs sta (t) = 7 p. 
yr 


(69) 


where, as for the Doppler case, Аза is a dimensionless 
parameter that can be related to the kinematic parameters of the 
transiting event (see Appendix C.6 for more details). In the 
dynamic limit, there is no simple parameterization of the 
Shapiro signal; therefore, we do not search for this signal. 
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Assuming a monochromatic PBH population, our goal is to 
derive a posterior distribution for the PBH mass fraction, 
вн = Орвн/Орм, as a function of the PBH mass, M: 
р(Љвн|ёѓ, M). We do this as follows: 


1. For each given value of вн and M, we use the Monte 
Carlo code developed by Lee et al. (20212) to generate a 
PBH population surrounding each of the pulsars in our 
array. From this distribution, we derive the amplitude of 
the static Doppler and Shapiro signals generated by the 
entire PBH population and the amplitude for the dynamic 
Doppler signal generated by the closest transiting PBH. 
Finally, we repeat this procedure for 2.5 x 10? realiza- 
tions to obtain the conditional probability distributions 
p(log,)Ar|fppy)» where J indexes pulsars in the array and 
A refers to any of the PBH signal amplitudes introduced 
in Equations (67), (68), and (69). In Figure 16 we report 
some of the distributions derived in this way.” 

. One at a time, we include the PBH signals given in 
Equations (67), (68), and (69) in the timing model, and 
we analyze our data to derive the posterior distributions 
for the various PBH signal amplitudes, p(log;, A|ót). 


21 From now on, we suppress the PBH mass, M, in the expressions for the 
conditional probabilities for the sake of notation. 
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Figure 15. The red solid (dashed) lines show our constraints at the 95% credible level for the two vector ULDM models considered in this analysis, for both the 
correlated (solid line) and uncorrelated (dashed line) regimes. The lower gray shaded regions correspond to regions of parameter space where the signal amplitude is 
smaller than the purely gravitational signal. The constraints from tests of the equivalence principle (“MICROSCOPE”) are shown in teal (Bergé et al. 2018), and the 
constraints previously set by the PPTA Collaboration are reported in yellow (Xue et al. 2022). 


Since the PBH signal in different pulsars is assumed to be 
independent, these distributions can be factorized as 


Np 
plogo Alt) = [| pogoA;lóD. 


1=1 


(70) 


Some of the р(108 4197) are reported in Figure 16. 
3. Finally, we can write 


Np 
DCfegulót) = П Ја logo Ar PCfegyullogio Ar) р (ор, А ді 
Т=1 


Np 
x [[ f41ogoA: ровіо Аг | /звн) P dog бі), 


1=1 


(71) 


where, in the second step, we used Bayes's theorem and 
assumed uniform priors оп log,, A; and вн. More details 
on each of these three steps can be found in 
Appendix C.6 or in Lee et al. (2021b). 


DM substructures can also possess macroscopic charges and 
interact with baryonic matter via long-range Yukawa interac- 
tions. These interactions can be modeled by a potential of the 
form 


Уша (0) = —à (72) 


GMMp e^, 
ГА 
where M and Mp are the masses of ће DM and pulsar, 
respectively, A is the range of the interaction, and & is the 
effective DM-baryon coupling, normalized against the 
gravitational coupling (also known as the Yukawa parameter). 
Here the DM can be either a particle or a macroscopic object 
such as a nugget of asymmetric DM (Detmold et al. 2014; Wise 
& Zhang 2014; Hardy et al. 2015; Krnjaic & Sigurdson 2015; 
Gresham et al. 2017, 2018). These Yukawa interactions can 
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Shapiro sta. | 


Doppler sta. | 
Doppler dyn У 


feen = 10? 


Figure 16. The black solid (dashed) lines show the posterior distributions 
P(log,yAstaldt) (p Mogi) Адуп|бі)) for a representative pulsar (J1909—3744). 
The filled distributions show the conditional probability distributions 
P (logo A|fpgu) for the same pulsar and different values of вн. In this plot, 
M = 107° (1079) Mo for the Doppler static (dynamic) signal and M = 10 Mo 
for the Shapiro static signal. 


arise from ап effective Lagrangian of the form 
L D gypXX + в qnn, where X and n are the effective DM 
and neutron fields, respectively, and @ is a massive (but 
potentially light) scalar or vector field. The effective coupling is 
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Figure 17. Constraints at the 9596 credible level on the local PBH abundance 
derived from the search for static Doppler (red shaded region) and static 
Shapiro signals (blue shaded region). The solid lines interpolate between the 
PBH masses simulated in this work, while the red dashed line shows an 
extrapolation of the contraints to higher masses. 


1076 1074 


Enx 
4тСтхти 
neutron mass. These interactions are constrained to be weaker 


than gravity for the mass range M < 100 GeV by the CMB, Lya 
forest (Xu et al. 2018), and direct detection experiments such as 
X-ray Quantum Calorimeter (XQC; Mahdawi & Farrar 2018) and 
Cryogenic Rare Event Search with Superconducting Thermo- 
meters (CRESST; Angloher et al. 2017) (for a review on these 
constraints see Xu & Farrar 2020). However, stronger-than- 
gravity fifth forces are allowed if М >> 100 GeV, even when X 
accounts for the entirety of the DM population. 

If present, these Yukawa interactions will contribute to the 
pulsar’s acceleration induced by a transiting DM substructure 
and contribute to the Doppler signal discussed before (the 
expression for the Yukawa contribution to the Doppler signal 
can be found in Appendix C.6). Therefore, as shown by 
Gresham et al. (2023), following a procedure similar to the one 
used to constrain the abundance of PBHs, we can constrain the 
value of the Yukawa parameter, & . Specifically, for each given 
value of & and M, we use the Monte Carlo code developed by 
Lee et al. (2021a) to generate a population of DM substructure 
surrounding each of the pulsars in our array. From this 
distribution, we derive the amplitude of the static Doppler 
signal generated by the closest transiting substructure by 
considering the acceleration induced by both the gravitational 
and Yukawa interaction. By repeating this procedure for 
multiple populations of DM substructure, we derive the 
distribution p(log,,A;|&). By plugging this quantity into an 
expression similar to the one given in Equation (71), we can 
derive p(&|ôt) and use this quantity to constrain бу. 


related to the coupling constants by & яз 


, Where т, is the 


6.2.2. Results and Discussion 


We start by searching for PBH signals on top of a GWB that 
we model as a power law with amplitude and spectral index 
allowed to float within the following prior ranges: 
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Figure 18. The 95% credible level for the fifth-force strength @ derived from 
the NGI5 data (red lines) is compared with constraints from NS kinetic heating 
(blue lines), equivalence principle constraints (green lines), and Bullet Cluster 
+ equivalence principle constraints (gray line). Solid (dashed) lines are 
deriving assuming А = 1 pc (А = 107! рс), while dashed-dotted (dotted) lines 
are derived assuming А = 10? pc (A= 10^? pc). 


0. 


logy Acws Є [-11, —18] and зсув Є [0, 7]. We find no 
statistically significant evidence for any of the PBH signals 
described in the previous section. Therefore, we proceed to set 
constraints on the local PBH abundance. The prior distributions 
used for the PBH signal parameters are reported in Table 3. 

The 95% upper limits on /рвн derived from the static 
Doppler and Shapiro signals are reported in Figure 17. The 
dynamic Doppler signal is too weak to produce any detectable 
signal for any of ће /рвн values considered. These are the first 
constraints on /рвн derived using real РТА data. As expected, 
our constraints are much weaker compared to the projections 
that were derived by Lee et al. (2021b) using mock data and 
including only white noise. Indeed, as already discussed by Lee 
et al. (2021b), the presence of a common red-noise process 
significantly reduces the sensitivity to PBH signals. 

Finally, in Figure 18 we show the constraints on à set by 
NGI25 data. These constraints are compared with several other 
constraints that can be placed оп ё. Specifically, in teal we 
show weak equivalence principle (WEP) constraints (Wagner 
et al. 2012; Shao et al. 2018; Sun et al. 2019; properly rescaled 
to take into account the finite range of the interaction, Gresham 
et al. 2023) derived by considering differential acceleration of 
baryonic test bodies toward the galactic center. In blue we 
report constraints from neutron star (NS) heating (assuming 
additional short-range DM-— baryon interaction) induced by 
DM capture (Gresham et al. 2023), derived from the coldest 
known NS to date—PSR J2144—3933 (Guillot et al. 2019). In 
gray we report the indirect constraints that can be derived by 
combining the fifth-force constraints on baryon—baryon 
interactions (Bergé et al. 2018; Fayet 2019) and Bullet Cluster 
constraints on DM—DM interactions (Spergel & Stein- 
hardt 2000; Kahlhoefer et al. 2014; see also Coskuner et al. 
2019; Gresham et al. 2023). We find that the NG15 constraints 
are competitive with WEP and NS kinetic heating, especially in 
the MZ 1072 Ms regime. However, the combined constraint 
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from the Bullet Cluster and MICROSCOPE dominates over all 
other constraints in the entire mass range of interest. We note, 
however, that the combined constraint is an indirect probe of 
the Yukawa parameter and depends on the specific form of the 
light mediator coupling to DM and baryons. Specifically, in 
deriving the combined constraints, we have assumed that the 
Yukawa DM-baryon interaction arises from a Lagrangian of 
the form £ D зх ФХХ + в onn. In addition, if only a 
subcomponent of DM interacts through the long-range fifth 
force, then the Bullet Cluster constraint will quickly lift, while 
the other constraints only deteriorate linearly with the DM 
fraction. 


7. Discussion 


The analysis of the NANOGrav 15 yr data set has produced 
the first convincing evidence for a stochastic background of 
GWs in the nanohertz frequency range. The origin of this 
background is still unknown. In this work, we considered 
various cosmological sources and compared them to the 
commonly studied astrophysical signal produced by a popula- 
tion of inspiraling SMBHBs. Specifically, we considered the 
signals produced by nonminimal inflation scenarios, SIGWs, 
cosmological phase transitions, several cosmic-string models, 
and domain walls. 

For each of these models, we identified regions of the 
parameter space that are compatible with the observed signal. 
We find that, with the exception of stable cosmic strings of field 
theory origin, all new-physics models considered in this 
analysis are capable of reproducing the GWB signal. Many 
models allow us in fact to fit the signal better than the SMBHB 
reference model, which is reflected in Bayes factors ranging 
from 10 to 100 (see Figure 2). When the new-physics signals 
are combined with the astrophysical one, we obtain comparable 
results. More precisely, in several models, the addition of the 
SMBHB signal leads to a slight decrease of the Bayes factor, 
which indicates that the SMBHB contribution does not help to 
improve the quality of the fit but merely increases the prior 
volume. In other models, on the other hand, such as PT-SOUND 
and DW-DR, adding the SMBHB signal on top of the new- 
physics signal can lead to a slight increase of the Bayes factor, 
indicating that the SMBHB signal can in fact play a dominant 
role in the total GW spectrum. For all four stable-string models, 
we find Bayes factors between 0.1 and 1. Cosmic superstrings, 
on the other hand, which are also stable but not of field theory 
origin, can explain the data at a level comparable to other new- 
physics sources. 

Despite the fact that some of the Bayes factors derived in this 
paper might suggest that a purely astrophysical interpretation of 
the signal is disfavored by the data, we caution against this 
interpretation. The Bayes factors do not account for the full 
range of uncertainties in both the cosmological and astro- 
physical signals and are prior dependent. It is conceivable that, 
as our understanding of SMBHB signals and our noise models 
improve, the tension between observations and astrophysical 
predictions will decrease, potentially weakening the evidence 
in favor of cosmological signals. 

Future data sets will improve the spectral characterization of 
the signal and improve our ability to discriminate cosmological 
sources from the SMBHB signal. Unfortunately, similarities in 
the spectral shape and theoretical uncertainties will make it 
challenging to definitively determine the origin of the back- 
ground using power spectrum characterization alone. However, 
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the observation of anisotropies could eventually resolve this 
debate, as the expected anisotropies generated by black hole 
binaries are significantly larger than those produced by most 
cosmological sources. Similarly, the detection of a continuous 
wave from a single binary would provide convincing evidence 
in favor of an astrophysical origin of the signal. Ultimately, 
measurements of the GWB spectral shape and angular power 
spectrum may be complemented by observations of its 
polarization content and possible deviations from Gaussian 
statistics, which can again help to discriminate between an 
astrophysical and a cosmological origin of the signal. 

It is worth emphasizing that in all parts of our analysis we 
described cosmological signals using effective parameters, e.g., 
Tx, аж, and H,.R,, for the phase transition models. Moving 
forward, it will be crucial to identify microscopic models that 
can reproduce the values of these parameters that we found to 
best fit the GWB signal. That is, in order to shed more light on 
the various cosmological interpretations of the signal, we 
require a better understanding of how the NANOGrav signal 
could possibly result from the fundamental parameters of a 
particle physics Lagrangian that describes the generation of 
GWSs in the early universe. 

Along with searching for a cosmological GWB, we also 
analyze our data to search for deterministic signals generated 
by models of ULDM and DM substructures. We do not find 
significant evidence for either of these signals. Nonetheless, we 
are able to place constraints on the parameters space of these 
models. For a wide range of ULDM models, our constraints 
compete or outperform laboratory constraints in the 
10 eV < т, < 10 ?eV mass window. The signal from 
DM substructures is harder to detect; as a consequence, we are 
able to set very weak constraints on the local abundance of 
these objects. Future data sets will improve our reach, but a 
better characterization of the GWB will be needed to probe 
realistic models of DM substructures. 

The discovery of a GWB will lead to significant break- 
throughs in our understanding of cosmology and particle 
physics. As future PTA data sets become available, we will 
establish the origin of the GWB. Regardless of whether the 
signal is of cosmological origin, we have shown how PTAs 
will undoubtedly contribute to exploring new physics, either as 
a discovery tool or as a new way to constrain the parameter 
space of BSM models. 
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Appendix A 
Parameter Ranges and Limits 


In this appendix, we specify the prior assumptions for all 
model parameters used in our analyses and report characteristic 
values for these parameters that we extract from the corresp- 
onding reconstructed 1D marginalized posterior distributions. 
In Table 2, we list our prior assumptions for the pulsar-intrinsic 
red-noise parameters Area and Yrea (Ag and Ya in Equation (4)), 
as well as for the SMBHB parameters Авнв and увнв. In the 
latter case, we work with a bivariate normal distribution for 
(logjyApup. Увнв) Whose mean and covariance matrix are 
given by 


NT 28  —0.026 
овнв = 10 ЗЕ 7 12 ) (AD) 


which we obtain by fitting the 102; Авнв and Эвнв distribu- 
tions extracted from the SMBHB simulations in the GWOn1y- 
Ext library in NGl5smbh (see Section 4). In Table 3 we list 
our prior assumptions for the model parameters of all new- 
physics models considered in this work, and in Table 4 we 
summarize various key values extracted from the corresp- 
onding 1D marginalized posterior distributions. Specifically, 
we report the Bayes estimator, including its standard deviation; 
the maximum posterior estimator, including the 68% credible 
interval; and (when applicable) the upper bound based on the 
requirement that the K ratio in Equation (10) should not be 
smaller than 1/10. 
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The Bayes estimator (0) of a parameter 0 with marginalized 
1D posterior probability distribution P(@|D, H) corresponds to 
the expectation value with respect to the distribution 
P(0/D,'H), while the standard deviation og of the Bayes 
estimator corresponds to the positive square root of the 
associated variance б, 


(0)— Ја 0 P(OID, H), о? = (82) — (0) 


=| f 4 8? P(0[D, m | = | / ав 0 P(OID, o. 
(A2) 


In practice, in a given analysis and for a given chain of MCMC 
samples, we compute the Bayes estimator and its standard 
deviation in terms of the corresponding sample mean and 
sample variance. The maximum posterior estimator Omax of a 
parameter 0 with marginalized 1D posterior probability 
distribution P(0|D, H) corresponds to the 0 value where 
P(6|D, H) reaches its global maximum across the predefined 
prior range, 


P(O@maxlD, H) = max P(6|D, H), (A3) 
9 
and the 68% Bayesian credible interval [Ө!", 0218 ] follows 
from integrating the posterior distribution P(6|D, H) over the 
regions of highest posterior density such that the integral 
returns an integrated probability of 68%, 


max 


І. арр, H) = 0.68, (A4) 
68 


where P(0|D, H) > Pg for all 0 є (om, дах) and some 
appropriate threshold Peg. Similarly, we can also construct 9596 
Bayesian credible intervals. Finally, we mention again that the 
K-ratio bound in the last column of Table 4 indicates the value 


Ox of the parameter 0 that returns К = E (see Section 3.2), 


көю = POR TO L L, (A5) 
РР] 6%, H) 10 
Note that, unlike all other quantities discussed in this section, 
the K ratio is not defined in terms of a posterior probability 
distribution, but rather in terms of a likelihood ratio, which 
makes it robust against our prior assumptions. 


Table 2 
Prior Distributions for the Pulsar-intrinsic Red-noise Parameters and the Parameters of the Astrophysical SMBHB Signal 


Prior Comments 


Pulsar-intrinsic Red Noise 


Parameter Description 
Ага Red-noise power-law amplitude 
^ed Red-noise power-law spectral index 


Log-uniform [—20, — 11] 
Uniform [0, 7] 


One parameter per pulsar 
One parameter per pulsar 


Supermassive Black Hole Binaries (SMBHB) 


(logi) Agus: Увнв) SMBHB signal amplitude and slope 


Могта(ивнв, ©внв) One parameter for PTA 


Note. The mean and covariance matrix of the Gaussian prior distribution for (log) Авнв; Увнв) аге given in Equation (А1). 
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Table 3 
Prior Distributions for the Parameters of the New-physics Models Considered in This Work 


Parameter Description Prior Comments 


Inflationary Gravitational Waves (IGW) 


Ta [GeV] Reheating temperature Log-uniform [—3, 3] One parameter per PTA 
r Tensor-to-scalar ratio Log-uniform [—40, 0] One parameter per PTA 
п, Tensor spectral index Uniform [0, 6] One parameter per PTA 


Scalar-induced Gravitational Waves (SIGW-DELTA) 


A Scalar amplitude Log-uniform [—3, 1] One parameter per PTA 
fa [Hz] Peak frequency Log-uniform [—11, —5] One parameter per PTA 
Scalar-induced Gravitational Waves (SIGW-GAUSS) 

A Scalar amplitude Log-uniform [—3, 1] One parameter per PTA 
fa [Hz] Peak frequency Log-uniform [—11, —5] One parameter per PTA 
А Width Uniform [0.1, 3] One parameter per PTA 
Scalar-induced Gravitational Waves (SIGW-BOX) 

A Scalar amplitude Log-uniform [—3, 1] One parameter per PTA 
Snin [Hz] Lower frequency Log-uniform [—11, —5] One parameter per PTA 
Јах [Hz] Upper frequency Log-uniform [—11, —5] One parameter per PTA 
Cosmological Phase Transition (PT) 

Т, [GeV] Transition temperature Log-uniform [—4, 4] One parameter per PTA 
Oy Transition strength Log-uniform [—2, 1] One parameter per PTA 
Н.К, Bubble separation Log-uniform [—3, 0] One parameter per PTA 
a Low-frequency slope (PT-BUBBLES) Uniform [1, 3] One parameter per PTA 

Low-frequency slope (PT-SOUND) Uniform [3, 5] One parameter per PTA 

b High-frequency slope (PT- Uniform [1, 3] One parameter per PTA 
BUBBLES) 

High-frequency slope (PT-SOUND) Uniform [2, 4] One parameter per PTA 

с Spectral shape width (PT-BUBBLES) Uniform [1, 3] One parameter per PTA 

Spectral shape width (PT-SOUND) Uniform [3, 5] One parameter per PTA 


Cosmic Strings (STABLE) 
Gu String tension Log-uniform [—14, —6] One parameter per PTA 


Metastable Cosmic Strings (META) 


Си String tension Log-uniform One parameter per PTA 
1-14, —1.5] 
JR Decay parameter Uniform [7, 9.5] One parameter per PTA 
Cosmic Superstrings (SUPER) 
Си String tension Log-uniform [—14, —6] One parameter per PTA 
P Intercommutation probability Log-uniform [—3, 0] One parameter per PTA 


Domain Walls (DW-SM) 


Т, [GeV] Transition temperature Log-uniform [—4, 4] One parameter per PTA 
Ax Energy fraction in domain walls Log-uniform [—3, 0] One parameter per PTA 
b High-frequency slope Uniform [0.5, 1] One parameter per PTA 
с Spectral shape width Uniform [0.3, 3] One parameter per PTA 
Domain Walls (DW-DS) 
Т, [GeV] Transition temperature Log-uniform [—4, 4] One parameter per PTA 
ANGE Amount of dark radiation Log-uniform One parameter per PTA 
[-3, —0.41] 

b High-frequency slope Uniform [0.5, 1] One parameter per PTA 
c Spectral shape width Uniform [0.3, 3] One parameter per PTA 
Ultralight Dark Matter (ULDM) 

Ai [s] ULDM signal amplitude Log-uniform [—9, —4] One parameter per PTA 
mg [eV] ULDM mass Log-uniform [—24, —19] One parameter per PTA 
e Earth normalized signal amplitude e" One parameter per PTA 

28 Pulsar normalized signal amplitude e" One parameter per 
pulsar* 

YE Earth signal phase Uniform [0, 27] One parameter per PTA 

Үр Pulsar signal phase Uniform [0, 27] One parameter per pulsar 


Primordial Black Holes (PBH-DYNAMIC) 


A Signal amplitude Log-uniform [—20, —12] One parameter per PTA 
To/ Tos Normalized time of closest Uniform [0, 1] One parameter per PTA 
approach 


Primordial Black Holes (PBH-STATIC) 
A Signal amplitude Log-uniform [—21, —13] One parameter per PTA 


Note. The asterisk indicates parameters that are present only in the uncorrelated ULDM analyses. 
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Table 4 
Bayesian Estimators, Maximum Posterior Values, and 6896 Credible Intervals for the Parameters of the New-physics Models 
Parameter Bayes Estimator Maximum Posterior 68% Credible Interval K Bound 
NP NP+SMBHB NP NP+SMBHB NP NP+SMBHB 
Inflationary Gravitational Waves (IGW) 
logi; Tin/GeV 0.02 + 1.60 —0.07 + 1.61 —0.53 —0.60 [- 1.51, 2.53] [- 1.89, 2.11] 
logor —14.06 + 5.82 —15.97 + 7.27 —10.14 —10.59 [-22.16, —6.58] [—23.03, —7.21] 
п, 2.61 + 0.85 2.68 + 0.97 2.02 2.08 [1.53, 3.92] [1.56, 4.03] 5.72 
logi; Авнв —15.60 = 0.56 —15.64 EE [— 16.20, —15.14] 
^/BHB 4.61 + 0.37 e 4.64 [4.26, 5.00] 
Scalar-induced Gravitational Waves (SIGW-DELTA) 
log, А —0.69 + 0.47 —0.71 + 0.49 —0.14 —0.17 [-— 1.00, —0.01] [—1.03, —0.02] 
log, of. /Hz —5.90 + 0.60 —5.93 + 0.67 -5 -5 [-6.17, —5*] [- 6.19, —5*] 
1089 ABHB —15.77 + 0.46 —15.71 [— 16.18, —15.29] 
^/BHB 4.65 + 0.35 EE 4.65 [4.31, 4.99] 
Scalar-induced Gravitational Waves (SIGW-GAUSS) 
1020 А —0.38 = 0.58 —0.36 = 0.61 —0.34 —0.29 [- 1.03, 0.20] [- 1.04, 0.24] 
logo f. /Hz —6.32 + 0.71 —6.30 + 0.73 —7.03 —6.85 [-7.25, —5.65] [—7.17, —5.57] 
А 1.35 + 0.70 1.30 + 0.70 1.60 1.54 [0.51, 2.07] [0.37, 1.92] 
102 0 ABHB — 15.72 + 0.46 — 15.65 e [— 16.14, —15.22] 
^/BHB 4.65 + 0.34 e 4.65 [4.32, 5.00] 
Scalar-induced Gravitational Waves (SIGW-BOX) 
log, А —1.06 + 0.52 —1.02 + 0.57 —1.26 —1.25 [—1.72, —0.82] [—1.68, —0.74] 
log Jimin /Hz —7.34 + 0.48 —7.29 + 0.55 —7.50 —7.48 [-8.01, —6.97] [—7.97, —6.84] 
logiofmax /Hz —6.06 + 0.65 —6.12 + 0.81 —5.40 —5.36 [—6.42, —5*] 1-6.45, —5*] 
logi; Авнв —15.70 + 0.49 — 15.64 [—16.15, — 15.21] 
^/BHB 4.65 + 0.35 EE 4.66 [4.31, 5.00] 
Cosmological Phase Transition (PT-BUBBLE) 
logi; T./GeV —0.76 = 0.49 —0.71 = 0.70 —0.90 —0.89 [—1.33, —0.39] [—1.34, —0.34] 
log, Qx —0.26 + 0.47 —0.23 + 0.52 Г 0.74 [0.03, #] [0.01, 1*] 
105 0 Н»Кх —0.42 + 0.26 —0.47 = 0.39 0 —0.06 [—0.56, 0*] [—0.58, 0%] 
a 2.04 + 0.48 2.07 + 0.49 1.97 2.01 [1.49, 2.54] [1.54, 2.63] 
b 1.97 + 0.58 1.98 + 0.58 г Г [I*, 2.32] [I*, 2.33] 
с 2.03 + 0.57 2.03 + 0.57 3 2.93 [1.69, 3*] [1.69, 3*] 
105 іо ABHB —15.68 + 0.51 —15.65 e [-16.17, —15.21] 
^/BHB 4.64 + 0.35 e 4.65 [4.30, 5.00] 
Cosmological Phase Transition (PT-SOUND) 
log Tx/GeV —1.84 + 0.41 —1.56 + 1.06 —2.00 — 1.95 [—2.33, — 1.48] [-2.31, — 1.30] 
log, Ox —0.22 + 0.44 0.14 + 0.56 —0.21 —0.15 [—0.37, 1*] [—0.34, 0.73] 
10210 HR. —0.81 = 0.36 —0.87 = 0.51 — 1.05 — 1.01 [— 1.28, —0.57] [—1.26, —0.45] 
a 3.58 + 0.47 3.74 + 0.54 3 3 [3*, 3.72] [3*, 3.98] 
b 2.87 + 0.57 2.92 + 0.57 2 2 [2*, 3.17] [2*, 3.25] 
c 4.16 + 0.56 4.09 + 0.57 5 5 [3.87, 5*] [3.77, 5*] 
logi; Авнв — 15.45 + 0.55 — 15.39 e [— 16.04, — 14.94] 
^/BHB 4.63 + 0.38 er 4.67 [4.27, 5.03] 
Cosmic Strings: Cusp-only Spectrum (STABLE-C) 
1020 Си —10.15 + 0.16 -1141 + 1.25 —10.18 —10.22 [—10.33, —10.01] [—12.13, —9.88] —9.67 
log, Авнв — 14.95 + 0.58 — 14.56 [—15.58, — 14.31 
^/BHB 4.34 + 0.38 e 4.24 [3.91, 4.66] 
Cosmic Strings: Kink-only Spectrum (STABLE-K) 
logo би —10.33 + 0.15 —11.34 + 1.17 —10.36 —10.38 [- 10.50, — 10.21] [-11.84, — 10.04] —9.87 
102 іо ABHB —15.04 = 0.61 — 14.56 [- 15.68, — 14.34 
^/BHB 4.39 + 0.38 e 4.28 [3.95, 4.72] 
Cosmic Strings: Monochromatic Emission (STABLE-M) 
logo би —10.53 + 0.14 — 11.47 = 1.09 —10.56 — 10.60 [— 10.68, — 10.42] [—11.91, — 10.27 —10.10 
102 іо ABHB —15.05 + 0.61 — 14.58 [— 15.67, — 14.34 
Увнв 4.39 + 0.38 e 4.28 [3.96, 4.73] 
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Table 4 
(Continued) 
Parameter Bayes Estimator Maximum Posterior 68% Credible Interval K Bound 
Cosmic Strings: Numerical Spectrum (STABLE-N) 
logo Си — 10.18 + 0.15 — 11.34 + 1.23 —10.21 —10.25 [—10.35, —10.05] [— 11.99, —9.90] —9.71 
1089 Авнв — 14.99 + 0.59 -14.56 [— 15.61, — 14.32] 
^/BHB 4.37 + 0.38 4.26 [3.93, 4.69] 
Metastable Cosmic Strings: Loops Only (METAL) 
logo би —5.80 + 0.78 —5.9 + 1.2 —5.04 —5.05 [—6.14, —4.84] [—6.19, —4.83] 
Ук 7.95 + 0.13 7.95 + 0.18 7.85 7.84 [7.81, 8.00] [7.80, 8.00] 
1089 Авнв —15.73 = 0.48 —15.67 et [— 16.18, —15.25] 
^/BHB 4.64 + 0.35 4.66 [4.30, 5.00] 
Metastable Cosmic Strings: Loops and Segments (META-LS) 
logo би —5.62 + 1.01 —5.70 = 1.40 —4.46 —4.44 [—6.26, —4.24] [—6.33, —4.15] 
VK 7.83 = 0.18 7.82 + 0.23 7.61 7.61 [7.59, 7.92] [7.57, 7.93] 
1089 Авнв — 15.71 = 0.49 —15.67 e [— 16.17, —15.23] 
^/BHB 4.64 + 0.35 4.65 [4.30, 5.00] 
Cosmic Superstrings (SUPER) 
logo би — 11.67 + 0.32 — 11.68 + 0.35 — 11.94 — 11.96 [— 12.08, — 11.50) [— 12.09, — 11.49] —9.97 
log; P —2.23 = 0.55 —2.21 + 0.57 -3 -3 [-3*, —2.01] [-3*, —1.99] 
108 о ABHB — 15.76 + 0.46 —15.70 [— 16.18, — 15.28] 
^/BHB 4.64 + 0.35 4.65 [4.30, 4.99] 
Domain Walls (DW-SM) 
logi; Тх/Сем —0.73 = 0.21 —0.65 = 0.49 —0.79 —0.78 [—0.96, —0.56] [-0.97, —0.51] 
log, Ож — 0.88 + 0.21 —0.87 = 0.32 —0.92 —0.90 [—1.10, —0.71] [—1.10, —0.66] 
b 0.74 + 0.14 0.74 + 0.14 0.5 0.5 [0.5*, 0.83] [0.5*, 0.83] 
с 2.01 + 0.70 1.92 + 0.74 3 3 [1.72, 3*] [1.57, 3*] 
1089 Авнв —15.60 + 0.50 —15.49 e [—16.06, —15.08] 
^/BHB 4.66 + 0.35 4.67 [4.32, 5.02] 
Domain Walls (DW-DR) 
logi; /GeV —0.98 + 0.15 —0.62 = 1.37 —0.94 —0.94 [—1.10, —0.82] [—1.27, —0.67] 
logo АЛ т —0.48 + 0.06 —0.87 + 0.71 —0.41 —0.41` [—0.49, —0.41*] [—0.63, 0.41*] 
b 0.74 + 0.14 0.74 + 0.14 0.5 0.5 [0.5*, 0.97] [0.5*, 0.83] 
с 1.95 + 0.68 1.83 + 0.73 3 3 [1.62, 3*] [1.44, 3*] 
1089 Авнв —15.33 + 0.65 —15.59 e [— 16.03, —14.40] 
^/BHB Ue 4.5] + 0.39 4.52 e [4.10, 4.90] 


Note. Values annotated with an asterisk are at the boundary of the prior range used in the analysis. 


Appendix B 
Median GW Power Spectra 


In Figures 3 and 4 in the main text, as well as in Figures 19 
and 20 in this appendix, we show median GW power spectra 
for all of the new-physics models considered in this work. The 
purpose of this appendix is to explain how these spectra are 
constructed. The main idea is to take the output of our Bayesian 
analysis, i.e., the reconstructed posterior distributions in the 
parameter space of the respective models, and map these 
distributions onto distributions of і ow values at each 
frequency in the NANOGrav frequency band. In practice, this 
means that we sample model parameter values from our 
posterior distributions and use these parameter values to 
evaluate the GW power spectrum at one frequency after 
another in small steps in the range from f= 1/7, to 
f=30/Tops. At each fixed frequency, we thus obtain an 
induced posterior distribution of fe se values, from which we 
can derive point values and uncertainty estimates. The median 
Осуу value of a distribution at fixed f defines the value of the 
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median GW spectrum at this frequency. Similarly, the equal- 
tailed 68% and 95% intervals around the median values 
provide an uncertainty estimate for the median GW spectrum; 
these intervals are shown in terms of the blue and red bands in 
Figures 19 and 20. 

We stress that the median GW spectra in Figures 3, 4, 19, 
and 20 represent effective power spectra that encapsulate global 
properties of the underlying distributions of actual GW spectra. 
In general, no individual GW spectrum at a given point in 
parameter space will exactly coincide with the median GW 
spectrum. For most models, the difference between the 
individual GW spectra and the median GW spectrum is rather 
mild. However, for some models, there can be noticeable 
differences, such as in the case of the SIGW models. Note, e.g., 
that the median GW spectrum of the SIGW-DELTA model in 
Figure 4 does not have the characteristic double-peak structure 
that each individual GW spectrum in this model has. The 
reason for this is clear: the median GW spectrum follows from 
a distribution of many individual GW spectra whose peaks are 
located at different frequencies. The double-peak structure of 
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Figure 19. Median GWB spectra (solid lines) for all new-physics models considered in this work except for the cosmic-string models, together with their 6896 and 
95% posterior envelopes. Median GWB spectra for the cosmic-string models are shown in Figure 20. In the left column (blue shading), we show the median GWB 
spectra for the new-physics models alone; in the right column (red shading), we combine the new-physics signals with the signal from SMBHBs. The gray violins are 
symmetrical representations of the 1D marginalized posterior probability density distributions of the GW energy density at each sampling frequency of the data. The 
dashed black lines show the GWB spectrum produced by inspiraling SMBHBs (see caption of Figure 3). 
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Figure 20. Same as Figure 19, but for the cosmic-string models considered in this work. See Appendix B for details. 


the GW signal in the SIGW-DELTA model is therefore washed 
out owing to the averaging over many individual GW spectra. 

Another caveat is that median GW spectra violating an 
experimental bound (e.g., the LVK bound) do not automati- 
cally indicate that the corresponding model is ruled out. Again, 
as the median GW spectrum is constructed from a distribution 
of individual GW spectra, the violation of an experimental 
bound merely signals that some (maybe most) individual GW 
spectra are in conflict with the experimental data. It is, 
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however, well possible that some fraction of the underlying 
distribution of individual GW spectra is not ruled out and is 
perfectly consistent with all experimental constraints. This is, 
e.g., true for the GW spectra from cosmic superstrings. The 
median GW spectrum of the SUPER model violates the LVK 
bound (see Figure 4). However, at the level of the model 
parameter space, this merely means that some parameter 
regions are experimentally ruled out, while other regions 
remain viable (see Figure 11). Another example is the median 
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GWB spectrum of the IGW model, which appears to violate the 
Моє bound if it is extended to high frequencies beyond the 
NANOGrav band (see Figure 4). However, the IGW model as a 
whole is not ruled out, as one can still choose the number of e- 
folds during reheating, М, for each individual GW spectrum. 
In some regions of parameter space, the М bound then 
persists to pose a problem, despite this additional parametric 
freedom and even in the limit № — 0; other regions, however, 
remain phenomenologically viable (see Appendix С.1 and 
Figure 22). 


Appendix С 
Supplementary Material 


С.1. Cosmic Inflation 


In Figure 22, we present constraints on the parameter space 
of IGWs, i.e., the IGW model discussed in Section 5.1. This 
parameter space is spanned by four quantities: the tensor-to- 
scalar ratio r, the spectral index of the primordial tensor power 
spectrum п, the reheating temperature Ти, and the number of e- 
folds during reheating М. However, because of the strong 
covariance between n, and г in Figure 5, we are able to 
eliminate r and work on the 3D hypersurface in parameter 
space where r is given by the linear fit in Equation (22). The 
effective parameter space spanned by Tin, п, and М is then 
subject to two constraints: (1) the upper limit on the allowed 
amount of extra relativistic radiation, ДМ, є at the time of BBN 
and CMB decoupling (see Equation (23)), and (ii) the upper 
limit on the amplitude of the SGWB set by the LVK 
Collaboration based on their first three observing runs (see 
Equation (24)). The № bound can be saturated for a critical 
value of the cutoff frequency fena in the GW spectrum, which, 
for given Tin, is solely controlled by the Hubble rate at the end 
of inflation, Heng. This means that ће Neg bound can be 
translated to a maximally allowed cutoff frequency ру’, which 
can then be converted to a maximally allowed Hubble rate 

«^ and ultimately an upper limit №“ on the allowed 
number of e-folds during reheating. Similarly, if the predicted 
GW signal at fj, = 25 Hz should exceed the ГУК bound, we 
can derive upper limits HI? and №" by requiring f.,4 not to 
exceed the lower end of "ihe LVK bid. for definiteness, we 
use fia —20 Hz. 

Our results for obtained in this analysis are presented 
in Figure 22 in terns of black and maroon contour lines, where 
the black solid lines refer to the №: bound, while the maroon 
dashed lines refer to the LVK bound. In the left panel we show 
our results for Np” itself, while in the right panel we compare 

пах to the diiration of reheating that one would naively expect 
in ит of standard single-field slow-roll inflation (see 


Equations (15) and (16)), 
end. 


2 тах A 
a =h] М = —In 


"E 


min 
end 


т? 1/2 
(Za) Мы, 


where H™ denotes the minimal Hubble rate that is necessary 
to realize the desired value of Tn after inflation, 


th \1/2 
vie y^ ті 
90 Mp 


H naive __ 


end 


(C1) 


min — 
end 


(C2) 
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In other words, for given Т, HZ? 
in the limit of instantaneous reheating, М, 


maximally allowed number of e-folds 


represents the Hubble rate 
th = 0. In terms of the 


m^ and the naive 


naive 


expectation Np , we can then construct 


= NANG = 2 end 
NU 


which is the quantity shown in the right panel of Figure 22. In 
addition to the contour lines for № and AN, we also 
display the samples in the 735-7; plane that we obtain from our 
MCMC chains. The density of these blue points approximates 
the 2D marginalized posterior probability density for the model 
parameters T, and n, 

In view of the results in Figure 22, we can draw several 
conclusions: (1) Both the №. bound and the ГУК bound are 
— relevant at large values of the spectral index. At 

Ta ^ 10 ? GeV, the number of e-folds during reheating 15 
only constrained as long as п, 2 2.5, while at Ти ~ 10? GeV, 
the bounds оп М, only appear at n; > 1. At f 7» fin, these п, 
values translate to a slope of the GW spectrum a in the range 
from a^ —1 to a~0.5, which is in the regime where we 
expect the approximate expression in Equation (25) to be valid 
(see also the discussion in Kuroyanagi et al. 2015). (ii) The 
ГУК bound only depends on Tin. This follows from the fact 
that it is derived from the requirement у’ = 20 Hz, which is 
independent of г and п, Typically, ће ГУК bound on Nj is 
much weaker than the А bound on Min; only in the parameter 
region where the Neg bound begins to disappear, Ny — oo, 
does Ше ГУК bound become competitive. (iii) The №. bound 
is particularly strong at large reheating temperatures and large 
values of the spectral index, where № can become even 
negative. This region, which we indicate by a dark-gray 
shading, is phenomenologically not viable and hence ruled out 
(see also the regions labeled Мг in Figure 5). In the excluded 
region, we also find AN, < 0, which indicates that reheating 
lasting for the naive number of e-folds Му." will definitely 
result in a violation of the Neg bound. (iv) Farther away from 
the excluded region, the upper limits on Nj become weak very 
fast, МЕ >> 1. In realistic models of inflation and reheating, 
where we typically expect Np ~ OC --- 10), these bounds 
should be easy to satisfy. We therefore conclude that most of 
the region shaded in light gray, as well as the entire white 
region in Figure 22, can host viable realizations of the IGW 
model. This includes, in particular, scenarios with a low 
reheating temperature, Ти, ~ 10 ^^ GeV, and a large spectral 
index of the primordial tensor power spectrum, п, ~ 3 --- 4 

In Figure 21, we show an extended version of the corner plot 
in Figure 5 including the SMBHB parameters Авнв and ^внв. 


АМ = 


(C3) 


C.2. Scalar-induced Gravitational Waves 


The generation of SIGWs relies on enhanced scalar 
perturbations on small scales, which can simultaneously lead 
to the production of PBHs. The parameter space of the SIGW 
models considered in this work is therefore subject to the 
constraint that the PBH mass density must not exceed the relic 
density of DM, Ірвн < 1 (see Equation (32)). This appendix 
explains how we evaluate this bound on the SIGW model 
parameters. We use the Press-Schechter formalism for 
spherical collapse (Press & Schechter 1974; Carr 1975) and 
mainly follow Inomata et al. (2017b) and Ando et al. (2018). 
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Figure 21. Same as Figure 5, but including the SMBHB parameters Авнв and увнв. The black dashed lines in the three lower right panels show the prior distributions 
for Авив and увнв, i.e., one 2D Gaussian and two 1D Gaussian distributions. 
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Figure 22. Constraints on the parameters of IGWs, i.e., the IGW model discussed in Section 5.1. The black solid lines refer to the Neg bound, while the maroon dashed 
lines refer to the LVK bound. The left panel shows the results for Му", while the right panel compares № to the duration of reheating that one would naively 


expect in models of standard single-field slow-roll inflation. The blue points correspond to the sampled points used to obtain the 2D marginalized posterior probability 
density for the model parameters Ти, and n, in Section 5.1. See Appendix C.1 for details. 
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The main quantity of interest is /рвн, which denotes the total follows Gaussian statistics, (M) is given by 
fraction of PBH DM, М 
dó б 
so f (е) 
fan = Гам). (C4) к Jro) 20500) 
M 2 
c (M) óc 
ES exp 5 ; (C7) 
with PBH mass function f given as (Ando et al. 2018) т б, 20°(M) 
1/4 12 where б„ is the critical density perturbation for the PBH 
fM) = a BM) (23) | 0.12 a) . formation, whose exact value depends on the shape of the 
1.6 x 107? Д gT) Ормі? JAM curvature power spectrum and satisfies 0.4 < 6. < 0.6 (Musco 
(C5) et al. 2005, 2021), and SM) is the variance of the coarse- 


Here M is the mass of the PBH, T is the temperature at the time 
of PBH formation, y is the ratio between the PBH mass and the 
horizon mass, and (M) is the PBH production rate. In this 
work, we only focus on PBHs produced during the radiation- 
dominated era. We choose у = 0.356, following Choptuik 
(1993), Koike et al. (1995), and Niemeyer & Jedamzik (1998), 
and assume that the mass of the PBH is proportional to the 
horizon mass when k — aH, 


M = y My = Mol 


ay 4.2 x 106 Mpc-! Y 
11075 | 


k 
(C6) 


We can then express the PBH mass as a function of the 
temperature using the standard temperature-wavenumber 
relation in ACDM. Assuming that the density perturbation 
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grained density perturbation, which in the radiation domination 
era is given by 


16 
2(k) = 
(К) і 


4 
[2 (£) w(2) a.» Pu. c9 
q \k k 


з r è а = 
Here W(x) is some window function smearing over k and 


[sin (4/3 к) — (q/43K) cos (q/V3k)] 
(q/ 43k 


is the transfer function for the radiation-dominated era. In our 
work, we choose 0, = 0.45 as a fiducial value (Ando et al. 
2018; Wang et al. 2019) and a Gaussian window func- 
tion W(x) = e-*72, 

Our analysis is subject to a few caveats. The correct value of 
the critical density б, depends оп the choice of the power 
spectrum and the window function(Ando et al. 2018; 
Musco 2019; Young 2019; Escrivà 2020; Escrivà et al. 2020; 


Tq, К) = 3 (С9) 
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Gow et al. 2021; Musco et al. 2021; Dandoy et al. 2023) and 
the nonlinear relation between density perturbations and 
density contrast (De Luca et al. 2019; Young et al. 2019; 
Escrivà 2020). In this work, we use a single value of 6, for 
three different power spectra for computational reasons. We 
also disregard any corrections from the QCD equation of state, 
which are expected to be important at the frequency range 
probed by РТА observations (Abe et al. 2021; Escrivà et al. 
2022a; Juan et al. 2022; Dandoy et al. 2023; Musco et al. 
2023). 

The statistics of the primordial scalar perturbations also 
strongly affect the abundance of PBHs since, in general, 
enhanced scalar perturbations come with different levels of 
non-Gaussianity. In this work, we focus on Gaussian 
primordial scalar perturbations. With non-Gaussianities, the 
same PBH abundance can be produced with scalar perturba- 
tions whose amplitudes are orders of magnitude smaller than 
those produced with Gaussian scalar perturbations. Therefore, 
primordial non-Gaussianity can dramatically affect the PBH 
bounds (Young & Byrnes 2013; Garcia-Bellido et al. 2017; 
Nakama et al. 2017; Pattison et al. 2017; Franciolini et al. 
2018; Cai et al. 2019; De Luca et al. 2019; Young et al. 2019; 
Iacconi & Mulryne 2023). In the parameter region where 
SIGWs manage to explain the NANOGrav signal, this would 
notably aggravate the PBH overproduction problem, which 
means that sizable non-Gaussianities could completely rule out 
the SIGW interpretation of the signal. In this work, we also 
neglect evolutionary effects on the PBH mass function, namely 
accretion and mergers, where accretion effects are expected to 
be small for sub-solar-mass PBHs (Ali-Haimoud & Kamion- 
kowski 2017; Raidal et al. 2019; De Luca et al. 2020; 
Vaskonen & Veermäe 2020). 

In summary, the final PBH DM fraction is quite sensitive to 
the assumptions discussed above. The bounds presented in the 
main text should be taken with a grain of salt owing to the 
uncertainties in the computation of /рвн: 

In Figures 23 and 24, we show extended versions of the 
corner plots in Figures 6 and 7 including the SMBHB 
parameters Авнв and ^внв. 


C.3. Cosmological Phase Transitions 


In the phase transition analysis discussed in the main text, we 
allow the low-frequency spectral index to float despite that 
causality predicts a — 3. We do this to capture possible double- 
peak spectral features with our simple power-law parameter- 
ization. However, it is not clear whether or not a strong and fast 
phase transition like the one needed to explain the observed 
signal could produce such a double-peak structure (Hindmarsh 
et al. 2017; Hindmarsh & Hijazi 2019). Therefore, in Figure 25 
we report the results of a phase transition analysis where we 
assume a —3. Figure 25 shows the reconstructed posterior 
distributions for the parameters o, Tą, and Н„К» of the PT- 
SOUND and PT-BUBBLE models, both for the case where the PT 
is assumed to be the only source of GWs (blue contours) and 
for the scenario where we consider the superposition of the PT 
and SMBHB signal (red contours). For the analyses where the 
SMBHB signal is included, we also report the posterior 
distributions for Авнв and увнв. For the PT-SOUND model we 
notice very minor differences compared to the analysis in the 
main text. However, for the PT-BUBBLE model we notice how 
the posterior for Ту is peaked to slightly smaller values for the 
reasons explained in the main text. 
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In Figures 26 and 27, we report the posterior distributions for 
all the parameters of the phase transition models, including the 
spectral shape parameters a, b, and c that were excluded from 
Figure 8. As expected for the PT-BUBBLE model, the low- 
frequency slope is peaked around a~2, which is the 
reconstructed slope of the GWB signal, while the posteriors 
for b and c are approximately flat. For the PT-SOUND model, 
the posterior for a is peaked around the lower limit of the prior 
range a = 3, and there is also а mild preference for larger 
values of the width parameter, as this would flatten the 
spectrum close to the peak. 


C.4. Cosmic Strings 


In our discussion of stable cosmic strings in the main text, 
we only present the reconstructed marginalized 1D posterior 
distributions for the dimensionless cosmic-string tension, Gu 
(see Figure 9). The four strings-plus-SMBHBs models 
STABLE-C+SMBHB, STABLE-K+SMBHB, STABLE-M+SMBHB, 
and STABLE-N+SMBHB, however, feature three model para- 
meters in total: Gu, Авнв, and ygup. In this appendix, we 
complement the discussion in Section 5.4 and show the corner 
plots for these parameters (see Figure 28). A notable feature of 
these corner plots is that the plateau region at small values of 
Си, which we had already observed in Figure 9, now also 
appears in the form of flat directions in the 2D posterior 
distributions for Gu and Авнв, as well as for Си and "ggg. 
Meanwhile, the 2D posterior distributions for Авнв and увнв 
represent distorted versions of the 2D Gaussian prior distribu- 
tion that we employ in our analysis, peaking at large Авнв and 
small “pup, where SMBHBs yield the dominant contribution to 
the signal. All four corner plots in Figure 28 are qualitatively 
identical and only display slight numerical differences. 

Next, let us turn to metastable strings. In the main text, we 
discuss the number density of closed string loops for the case of 
stable strings (see Equation (43)), as well as for metastable 
strings (see Equation (46)). In the META-LS model, we need in 
addition the number density of string segments that form when 
long strings and closed loops break apart as a consequence of 
monopole nucleation. In order to compute the GW signal from 
segments, we use again Equations (40) and (42), where we 
replace the loop number density п, by the segment number 
density n, and the GW power spectrum P, in Equation (41) by 
P, —4/k, which was derived by Martin and Vilenkin in the 
approximation of a straight string segment in Martin & 
Vilenkin (1997). Meanwhile, all relevant expressions for the 
segment number density n, were computed in Leblond et al. 
(2009) and Buchmuller et al. (2021), which we shall 
summarize in this appendix. For more details, we refer 
to Leblond et al. (2009) and Buchmuller et al. (2021). 

First, we consider segments that form when long strings 
break apart. We denote their number density by A,, which we 
decompose into three different contributions that are relevant at 
different times and in different parameter regimes, 

fis(€, t) = nj (C, t) + + ANCE, t). (C10) 
Here л (£, t) is the number density of segments that originate 
from long strings, assuming the scaling regime to end during 
radiation domination, evaluated during radiation domination 
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Figure 24. Same as Figure 7, but including the SMBHB parameters Авнв and Эвнв- 
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Figure 25. Same as in Figure 5, but for the PT-SOUND (left panel) and PT-BUBBLE models with а low-frequency slope fixed to the value predicted by causality, 


ie., а= 3. 


and after the onset of the decay regime, 


I? (+1)? 


"dee vim 2 t) 


«| пее | 


where Лу (0,1) is the number density of segments that 
originate from long strings, assuming that the scaling regime 
ends during radiation domination, evaluated during matter 
domination and after the onset of the decay regime, 


gi 


ЕС +t) + 5 Tang — t)(t +3 2] 


t)O(t — t) 


exp 


i) 4 trone- +319), (Сп) 


як d. D) = e(t = teg) O (faq — ts) 


(C12) 


and Лу (£, t) is the number density of segments that originate 
from long strings, assuming that the scaling regime ends during 
matter domination, evaluated during matter domination and 
after the onset of the decay regime, 


T2 
teg) -4 exp 


m 


x EC t4 ; TONG — t) (t + »)] 


APP, т) = O(r — KOH — 


(C13) 


In view of these expressions, several comments are in order: (1) 
Throughout our analysis, we assume that GW emission by 
segments is as efficient as GW emission by loops, i.e., we work 
with Г=50 for both loops and segments. (ii) All three 
expressions depend on the dimensionless correlation length of 
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the long-string network, for which we use the attractor values 
in the VOS model, &, = 0.27 and Є, = 0.63 during radiation 
and matter domination, respectively. (iii) For a given choice of 
parameter values, йу, й, and п,“ never contribute 
simultaneously to the segment number density. For teq > f, 
only A; and Лу yield nonvanishing contributions (first й, at 
t< feq and then л, at t > teq), whereas for teq < t, the number 
density of segments from long strings is solely determined by 
по" at all times t> t, 

In addition to segments from long strings, there is also a 
population of segments that form when closed string loops 
begin to break apart because of monopole nucleation. We 
denote the number density of this population by i, 1, where the 
index refers to the fact that Ai, only describes the first 
generation of segments that form when closed loops break apart 
for the first time. This first generation then gives rise to a 
second generation of segments that follow from monopole 
nucleation on first-generation segments, and so on and so forth. 
We comment on these higher generations further below. Before 
that, however, we discuss A, 1, which we decompose again into 
three contributions, 


fii, t) = ANG, t) + ATE, т) + fg (t). (C14) 


Here ft) is the number density of first-generation 
segments that originate from closed loops that formed during 
radiation domination, evaluated during radiation domination 
and after the onset of the decay regime, 


fij (5, t) = Ot — ts) O(teqg — 1) O(teq — t) Га 


1 
х K - t) + = Гош – jore. 2), (C15) 
where й, 1 (2, t) is the number density of first-generation 
segments that originate from closed loops that formed during 
radiation domination, evaluated during matter domination and 


THE ASTROPHYSICAL JOURNAL LETTERS, 951:L11 (56pp), 2023 July 1 Afzal et al. 


EN PT-BUBBLE + SMBHB 
| PT-BUBBLE 


a 
= N N 
oo c 


b 
io B, oe 
пом 


.] 
a 


c 
к N N 
л Oo c 


10510 Авнв 
|. II] 
к e н 
чо Ct 


YBHB 


—2 0 —1 0 —2 —1 1.5 2.5 L5 25 15 25 -17 —15 4 5 
19510 Tx /GeV 19510 Q x 10810 Ax Rx a b c log10 ABHB YBHB 


Figure 26. Same as Figure 8, but including the spectral shape parameters a, b, c and SMBHB parameters Авнв and Увнв: 


after the onset of the decay regime, after the onset of the decay regime, 


ANG, t) = (t — 5)8(t& — t) Olt — teq) 


пт, f) = Ot — 5)8 — teq) O(teq — tx) Та Dy e азу Тена? 
x (C17) 
x K 5) 4 ; Gn ы pee. t), (C16) PE + TGp 0)? 
1 2 
and йуу (2,1) is the number density of first-generation {oar ea ь) + 2 ТОШЕ | 
segments that originate from closed loops that formed during + 0.45€ + Гби 0179? [E (r) — RW — Ft) + FI) 
matter domination, evaluated during matter domination and | і (C18) 
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Figure 27. Same as Figure 26, but for the PT-SOUND model. 


where F and F; are shorthand notations for the following 
expressions involving the hypergeometric function 5F;: 


Гав 
БО = oF ‚ -Вт+1 ; 
(x) => б B, —8;n B Tu 
п-1 n—f 
"DEC. *___ п=1,2, B=031. (C19) 
РАГСі n — 8 


These results were derived in Buchmuller et al. (2021) based on 
the loop number densities in the so-called BOS model (Blanco- 
Pillado et al. 2014), which explains the occurrence of factors 
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such as 0.27, 0.45, and 0.31 in Equations (C17) and (C18) (see 
Buchmuller et al. 2021 for more details). 

Finally, the expressions in Equations (C15), (C16), and 
(C17) enable us to estimate the total number densities of all 
generations of segments that descend from closed loops 
breaking apart because of monopole nucleation. In principle, 
these number densities are described by a partial integro- 
differential equation where the abundance of the nth segment 
generation acts as a source term for the (n + 1)th generation. 
Formally, this partial integro-differential equation can be 
solved analytically in terms of ап infinite recursive 
series (Buchmuller et al. 2021). The numerical evaluation of 
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Figure 28. Posterior distributions for the STABLE+SMBHB strings models including the SMBHB parameters Авнв and увнв. We also report the marginalized logi; Gu 


posterior distributions for the STABLE string models (blue lines). 


this series is, however, technically demanding, which is why 
we choose to follow a different approach for the purposes of 
this paper. As shown in Buchmuller et al. (2021), it turns out 
that the total number densities for segments from closed loops 
result in predictions for the GW spectrum that are very similar 
to those obtained from the corresponding first-generation 
number densities in Equations (C15), (C16), and (C17)—up 
to a numerical fudge factor of O(10). At the level of the number 
densities, it therefore suffices to rescale all first-generation 
number densities by a constant factor in order to obtain 
effective number densities for all generations of segments that 
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descend from closed loops breaking apart, 


й, (£, t) — fudge x fij, t), 

й, (£L, t) — fudge KANA, t), 
ñ (£, t) > fudge x йү (4, t). (C20) 
We stress that the functional dependence on / and t is typically 
not the same for the first-generation and total number densities; 
the rescaling in Equation (C20), however, achieves comparable 
results at the level of the GW spectrum. In our analysis, we 
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Figure 29. Same as Figure 10, but including the SMBHB parameters Авнв and Эвнв- 
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consistently use a fudge factor of 5, which is a characteristic 
value across large regions of parameter space. 

In Figures 29 and 30, we show extended versions of the 
corner plots in Figures 10 and 11 including the SMBHB 
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signals as 
GM м 2 я . =| ^ 
hp(t) — Y а. (41 + x5bp — sinh (x)) (C21) 
hs(t) = 2GM log(1 + x2), (C22) 


parameters Авнв and ^внв. 


C.5. Domain Walls 


In Figure 31, we report the posterior distribution for all the 
parameters of the domain wall models, including the spectral 
shape parameters b and c that were excluded from Figure 12. 
For both the DW-SM and DW-DR models we find a flat posterior 
for the high-frequency slope, b. This result is expected since 
most of the low frequency bins are fit by the low-frequency tail 
of the spectrum. For the width parameter, c, instead, both 
models prefer values near the upper range of our prior range. 
This preference for wider spectra derives from a mismatch 
between the reconstructed slope of the GWB, a — 2, and the 
one predicted by causality, а ~ 3. Larger width parameters 
make the spectrum from domain walls flatter near the peak and 
allow for a better fit to the data. 


C.6. Dark Matter Substructures 


In this appendix we provide more details on the procedure 
that we followed to derive the constraints on fpgy reported in 
Section 6.2. 

Given a PBH with position relative to the pulsar given by 
r(t) = го + vt, where ro and v are the initial PBH position and 
velocity, respectively, we can write the Doppler and Shapiro 


where xp = (t — tp o)/Tp and xs = (t — tso)/rs. These expres- 
sions only include cubic terms in time ¢ and have been 
previously derived in Dror et al. (2019) and Lee et al. (2021b). 
For the static limit in which т >> Tops, these expressions can be 
expanded in series of ю/т. Since ће O(t?) terms would be 
degenerate with the timing model, the measurable new-physics 
signal can then be parameterized as a term xf as 


A sta 
h(t)p(s, — Е (С23) 


where Арсѕ), ха for both the Doppler and Shapiro static signal 
cases are dimensionless amplitudes given by 


26M 5 1,0 1 z 
Ap, stat = уг d . | bp 


ard 


942, J 72 
| 1 1 28.0 / D s] (C24) 


3ту (l + io / ть)??? 


,4GM tso 3 — 180/75 
r 4 2 [23^ 
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For the Doppler case, in the dynamic limit when т < То, 
the dominant contribution would come from the first term in 


Equation (C21) where \/1 + хо œ |t — tol. Up to linear order 


As stat = (C25) 
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Figure 31. Same as Figure 12, but including the spectral shape parameters b, с and the SMBHB parameters Авнв and 7внв: 
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in xp, we can write 


hp,aya (t) = Ap, dyn(t = to)O(t - to), (C26) 
where Ap ауп is the dimensionless amplitude given by 
2GM > > 
Ap,dyn = E а. bp. (C27) 
| ту? 


Given the expressions in Equations (C21) and (C22) for the 
Doppler and Shapiro signals, we use the MC developed in Lee 
et al. (2021a) to derive the distributions p(log,,A/|fppy)- 
Specifically, we proceed as follows: 


1. For each pulsar we generate a population of Npgy PBHs, 
where Npp is implicitly defined by the relation 
ррмУ 


М№вн = вн , 
Мрвн 


(С28) 


where the simulation volume, V, is a sphere of radius 
К = ?Тъь» centered around the pulsar position for the 
Doppler signal and a cylinder with the same radius and 
height given by the Earth— pulsar distance for the Shapiro 
signal. Here v = 340 km s7! is the average PBH velocity, 
and Тоь is the observation time of the /th pulsar. 

. For each of these PBHs we generate a random initial 
position and velocity. Since PTA searches are only 
sensitive to DM subhalos in the neighborhood of the solar 
system, we expect the position distribution to be uniform. 
Therefore, we use the probability density function 
f(r) = 1/V to sample initial positions. To sample PBHs' 
velocity, we use a Maxwell-Boltzmann distribution with 
Vo = 325 km A Vese = 600 km s and the angular 
dependence assumed to be isotropic. 


N 


3. The PBHs’ signals are then classified as dynamic if they 
satisfy the condition Ту 7719» т and static 
otherwise. 


4. To evaluate Asta, we sum the static signals of all PBHs 
computed by using Equations (C21) and (C22), and we fit 
the resulting signal to a cubic polynomial in time to 
extract the г term. To compute Ap ауп, we take the closest 
transiting object and compute the signal amplitudes using 
Equation (C27). 


All the previous points are repeated for 2.5 x 10? realiza- 
tions to obtain the distributions р(108 о A;|fpgg). Given the 
conditional distributions p(logjA;|fpgg) and the posterior 
distribution p(log,,A7|6t) derived by analyzing our data, we 
can write 


Np 
рані) П Госел Arp ogy, Arló04 logo Ar. 
1=1 
(С29) 
Then, using Bayes's theorem, we can rewrite 


(C30) 
p (logo Аг) 


РО(Рвн Овід І) = 


Our priors on p(fppy) and p(log,,A;) are uniform, which 
allows us to eventually write Equation (C29) as 


Np 
P (УРвніді) сс П Jp ogi Ar] enu) 
І-І 


pogoA;jlót) d log, А, (С31) 
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where the œ implies that the p(fpgy|dt) would be subject to the 
normalization condition, Í, T P(fppylOt) вн = 1. 

In the presence of а DM— baryon fifth force in the form of a 
Yukawa potential in Equation (72), assuming point-mass DM, 
the pulsar frequency shift due to Doppler effect is given by 
(Gresham et al. 2023) 


(à) - SGM 7 1 
ыы аа t f 213/2 
V / fifth тру (1 + xp) 


x (1 af zm ES xB )е eri GS + xp¥)dxp. 


The integral in Equation (C32) has to be computed numeri- 
cally, and the signal due to the fifth force can be computed by 
performing an additional integration over time and subtracting 
away degeneracies with timing model parameters. The total 
signal is the sum of the gravitational and the fifth-force 
contribution, Ap totalt) = Ap sen) + Йразу(). In this analysis of 
the fifth-force constraint, we only consider the scenario where 
the DM substructure makes up the entirety of the DM local 
density, which is equivalent to taking јьвн = 1 for the 
gravitational contribution. Parameterizing the signal as 
Рр total) = ‚мар similar to the PBH case, the amplitude 
Арюош 15 a random variable dependent on A and а. The 
probability distribution function Р(109 Ар, юыА, Œ) can Бе 
extracted again by Monte Carlo simulations and Bayes's 
theorem. Finally, the posterior distribution of & and A, 
P(à, A|ót), is given by the analog of Equation (C31), 


(C32) 


Np 
p(&, М0) ox П [paos Ала, A)p(log,, A;|ót) d 108,4. 
1=1 


(C33) 
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